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Introduction 

The  ‘Time-resolved  Spectral  Optical  Breast  Tomography”  research  project  aims  to  develop  a  near  real¬ 
time  three-dimensional  (3D)  spectral  tomographic  imaging  algorithm  with  use  of  the  cumulant 
approximation  to  radiative  transfer  to  model  light  propagation  in  tissues.  This  project  in  this  reporting 
period  involves  theoretical  modeling  of  photon  migration  in  tissues  and  image  reconstruction,  and 
working  with  the  experimental  group  to  apply  and  test  the  image  reconstruction  algorithm  using 
experimental  data  obtained  from  phantoms.  Significant  advances  were  made  during  the  current  reporting 
period. 


Body 

The  tasks  performed  and  the  progress  made  during  the  current  reporting  period  include  theoretical 
modeling  of  photon  migration  in  tissues  and  image  reconstruction,  and  working  with  the  experimental 
group  to  apply  and  test  the  image  reconstruction  algorithm  using  experimental  data  obtained  from 
phantoms. 

Theoretical  modeling  of  photon  migration  in  tissue  and  image  reconstruction 

We  continued  to  enhance  the  3D  tomographic  image  reconstruction  algorithm  based  on  the  new  cumulant 
transport  model[l]  to  make  use  of  spectral  information  when  observations  of  multiple  wavelengths  are 
available  (Task  1 .3).  By  scanning  a  point  source  on  the  grids  of  the  input  plane  of  a  slab  and  measuring 
light  intensities  on  a  detector  array  on  the  exit  plane  of  the  slab,  a  set  of  four-dimensional  (4D)  data  is 
formed. [2,  Appendix  1]  The  spectral  information  adds  an  additional  dimension  of  the  data.  We  started  to 
investigate  the  optimal  approach  to  analyze  this  huge  dataset  (Task  1 .5).  Some  preliminary  results  were 
obtained  from  the  application  of  information  theory  to  the  simulated  and  experimental  dataset  using 
Independent  Component  Analysis  (ICA).[3,  4,  Appendix  4]  Improvement  in  the  quality  of  image 
reconstruction  was  observed. 

We  improved  the  3D  tomographic  image  reconstruction  algorithm  to  use  a  L-curve  method  guided  by  the 
signal-to-noise  ratio  of  the  dataset  to  determine  the  regularization  parameter  (Task  1 .4). 

We  also  studied  the  nonlinear  effect  of  the  multiple  passages  of  light  through  an  absorption 
inhomogeneity  for  optical  imaging  (Task  1).  We  derived  the  nonlinear  correction  factor  (NCF)  using  the 
cumulant  solution  to  radiative  transfer.  NCF  was  verified  and  supported  by  both  Monte  Carlo  simulations 
and  experiments.  The  nonlinear  correction  using  NCF  was  shown  to  correct  the  underestimation  of 
absorption  by  the  conventional  linear  perturbation  scheme  of  optical  imaging  when  the  inhomogeneity  is 
strong  in  absorption.[5,  6,  Appendix  3,5] 

Apply  and  test  the  image  reconstruction  algorithm  using  experimental  data 

We  worked  with  the  experimental  group  at  the  Institute  of  Ultrafast  Spectroscopy  and  Lasers  at  the  City 
College  of  New  York.  The  image  reconstruction  algorithm  was  modified  to  accept  the  data  from  the 
experimental  group  (Task  2).  Experiments  were  performed  to  image  objects  inside  tissue-like  Intralipid- 
10%  suspensions  in  water.  Experimental  results  were  analyzed  and  images  were  reconstructed  using  the 
3D  tomographic  image  reconstruction  algorithm  (Task  2.2).[3,  4,  Appendix  4]  We  are  working  to  fine 
tune  the  algorithm  using  the  experimental  data. 

Further  experiments  will  be  performed  on  breast  phantoms.  Experiments  to  use  multiple  wavelengths  and 
different  configurations  (transmission  and  backscattering)  of  the  setup  will  be  used  to  test  and  improve 
the  3D  tomographic  image  reconstruction  algorithm. 
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Key  Research  Accomplishments 

•  Developed  and  enhanced  the  3D  tomographic  image  reconstruction  algorithm  to  use  a  L-curve 
method  guided  by  the  signal-to-noise  ratio  of  the  dataset  to  determine  the  regularization 
parameter. 

•  Developed  and  enhanced  the  3D  tomographic  image  reconstruction  algorithm  to  make  use  of 
spectral  information  in  multiple  wavelength  measurements. 

•  Developed  a  novel  information  theory  approach  to  analyze  the  dataset  for  image  reconstruction. 

•  Extended  the  range  of  applicability  of  the  linear  inversion  scheme  for  optical  imaging  using  a 
nonlinear  correction  factor  for  inhomogeneities  strong  in  absorption 


Reportable  Outcomes 

Journal  Papers: 

1 .  Cai,  W.,  M.  Xu,  and  R.R.  Alfano,  Three  dimensional  radiative  transfer  tomography  for  turbid 
media.  IEEE  JSTQE,  2003.  9:  p.  189-198.  (Appendix  1) 

2.  Xu,  M.  and  R.R.  Alfano,  More  on  patterns  in  Mie  scattering.  Opt.  Comm.,  2003.  226(1-6):  p.  1- 

5.  (Appendix  2) 

3.  Xu,  M.,  Light  extinction  and  absorption  by  arbitrarily  oriented finite  circular  cylinders  using 
geometrical  path  statistics  of  rays.  App.  Opt.,  2003.  42:  p"  6710-6723 

4.  Xu,  M.,  W.  Cai,  and  R.R.  Alfano,  Multiple  passages  of  light  through  an  absorption 
inhomogeneity  in  optical  imaging  of  turbid  media.  Opt.  Lett.,  2004  (in  press)  (Appendix  5) 

Presentations  and  Proceeding  Papers: 

5.  Xu,  M.,  W.  Cai,  and  R.R.  Alfano.  Nonlinear  multiple  passage  effects  on  optical  imaging  of  an 
absorption  inhomogeneity  in  turbid  media,  in  European  Conference  on  Biomedical  Optics: 
Photon  migration  and  Diffuse-light  imaging.  2003.(Appendix  3) 

6.  Xu,  M.,  et  al.,  Simulated  and  experimental  separation  and  characterization  of  absorptive 
inhomogeneities  embedded  in  turbid  media,  in  Biomedical  topical  meetings  on  cd-rom  (osa). 
2004:  Fontainebleau  Hilton  Resort  and  Towers,  Miami  Beach,  Florida.  WF25  (Appendix  4). 

7.  Al-rubaiee,  M.,  et  al.,  Time-resolved  and  quasi-continuous  wave  three-dimensional  tomographic 
imaging,  in  Femtosecond  laser  applications  in  biology.  2004:  Palais  de  la  Musique  et  des 
Congres  de  Strasbourg,  Strasbourg,  France. 

Grant  application: 

Applied  for  breast  cancer  concept  award  “Localization  and  Identification  of  Tumor  for  Optical 
Breast  Imaging”  for  BC03-CA. 


Conclusions 

The  work  carried  out  during  the  current  reporting  period  builds  on  and  affirms  some  of  our  earlier 
inferences  and  leads  to  the  following  conclusions.  First,  the  L-curve  method  guided  by  the  signal-to-noise 
ratio  of  the  dataset  and  the  information  theory  approach  using  independent  component  analysis  to  analyze 
the  dataset  were  found  to  improve  the  quality  of  image  reconstruction.  Second,  the  3D  tomographic  image 
reconstruction  algorithm  was  enhanced  to  make  use  of  spectral  information  in  multiple  wavelength 
measurements;  Third,  the  correction  for  the  nonlinear  effect  of  the  multiple  passages  of  an  absorption  site 
by  light  was  shown  to  be  essential  in  optical  imaging  to  characterize  properly  inhomogeneities  strong  in 
absorption.  Fourth,  the  theoretical  formalism  and  computer  algorithm  for  3D  tomographic  image 
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reconstruction  shows  (with  simulated  and  experimental  data)  the  potential  to  provide  fast  3D  images  of 
the  scattering  and  absorption  objects  at  various  depths  in  turbid  media. 
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Three-Dimensional  Radiative  Transfer 
Tomography  for  Turbid  Media 

W,  Cai,  M.  Xu,  and  R.  R.  Alfano 


Abstract — The  photon  distribution,  as  a  function  of  position, 
angle,  and  time,  is  computed  using  the  analytical  cumulant  solution 
of  the  Boltzmann  radiative  transfer  equation  (RTE).  A  linear  for¬ 
ward  model  for  light  propagation  in  turbid  media  for  three-dimen¬ 
sional  (3-D)  optical  tomography  is  formed  based  on  this  solution. 
The  model  can  be  used  with  time  resolved,  continuous  w  ave  (CW), 
and  frequency-domain  measurements  in  parallel  geometries.  This 
cumulant  forward  model  (CFM)  is  more  accurate  than  that  based 
on  the  diffusion  approximation  of  RTE.  An  inverse  algorithm  that 
incorporates  this  CFM  is  developed,  based  on  a  fast  3-D  hybrid- 
dual-Fourier  tomographic  approach  using  multiple  detectors  and 
multiple  sources  in  parallel  geometries.  The  inverse  algorithm  can 
produce  a  3-D  image  of  a  turbid  medium  with  more  than  20  000 
voxels  in  1-2  min  using  a  personal  computer.  A  3-D  image  recon¬ 
structed  from  simulated  data  is  presented. 

Index  Terms — Absorption  and  scattering,  forward  model,  in¬ 
verse  algorithm,  optical  tomography,  photon  migration,  radiative 
transfer  equation  (RTE). 

I.  Introduction 

OVER  THE  PAST  decade,  optical  tomography  has  been 
investigated  as  a  noninvasive  imaging  method  that  uses 
nonionizing  near-infrared  (NIR)  light  to  obtain  images  of  the 
interior  of  the  breast.  Unlike  X-ray,  which  is  attenuated  through 
media  by  ionizing  the  electrons  at  inner-orbits  of  atoms,  NIR 
light  uses  the  vibrational  overtones  for  different  molecular  com¬ 
ponents  in  the  structures  of  tumor.  NIR  light  may  be  used  to 
create  image  based  on  the  molecular  change,  which  may  be  used 
to  improve  sensitivity  and  specificity  in  the  early  diagnostics 
of  breast  cancer.  Breast  tissues  scatter  light  strongly,  and  blur 
the  direct  shadow  image  of  a  tumor.  A  technique,  known  as  in¬ 
verse  image  reconstruction,  has  been  investigated  to  overcome 
the  problem  of  multiple  scattering.  Some  obstacles  in  the  de¬ 
velopment  of  optical  tomography  are  inaccuracy  of  the  com¬ 
monly  used  diffusion  forward  model,  and  lack  of  a  fast  inverse 
algorithm  able  to  realize  a  three-dimensional  (3-D)  image  re¬ 
construction  of  a  breast  for  clinical  use. 

One  critical  issue  is  the  forward  model,  which  should  cor¬ 
rectly  simulate  photon  propagation  in  the  medium.  The  most 
commonly  used  forward  models  were  built  based  on  solution 

Manuscript  received  November  6,  2002;  revised  February  11,  2003.  This 
work  was  supported  in  part  by  the  U.S.  Army  Medical  Research  and  Materiels 
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the  Office  of  Naval  Research  (ONR). 

The  authors  arc  with  the  Institute  for  Ultrafast  Spectroscopy  and  Lasers,  New 
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and  Applications,  Department  of  Physics,  The  City  College  and  Graduate 
Center  of  City  University  of  New  York,  New  York,  NY  10031  USA  (e-mail: 
alfano@scisun.sci.ccny.cuny.edu). 
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of  the  diffusion  equation,  which  is  the  lowest  approximation 
of  the  radiative  transfer  equation  (RTE)  [l]-[5].  The  forward 
models  based  on  the  diffusion  approximation  (DA)  give  a  large 
error  when  the  distance  d  between  a  voxel  and  a  source  is 
small.  Furthermore,  the  photon  distribution  still  maintains  a 
strong  anisotropy  in  a  deeper  region  away  from  a  source,  which 
will  be  shown  later  in  this  paper.  Unfortunately,  contributions 
from  near  surface  voxels  to  measured  signals  are  often  larger 
than  contributions  from  the  voxels  deep  inside  the  medium. 
Inaccuracy  of  the  DA-based  forward  model  may  lead  to  a 
failure  in  image  reconstruction,  especially  for  small  hidden 
objects  deep  inside  the  medium.  The  total  weight  matrix  should 
be  inverted.  The  large  elements  in  the  matrix,  which  play  a 
more  important  role  in  inversion,  are  evaluated  incorrectly  in 
DA  models.  The  shortcoming  of  DA  is  well  recognized,  but  it 
is  still  broadly  applied  due  to  the  difficulty  in  directly  solving 
the  radiative  transfer  equation.  Hielscher  et  al  [6]  and  Vihunen 
et  al.  [7]  developed  numerical  solutions  of  RTE  for  optical 
tomography. 

Recently,  we  have  developed  an  analytical  solution  of 
RTE,  based  on  cumulant  expansion,  in  an  infinite  uniform 
medium  with  an  arbitrary  phase  function  [8],  [9].  It  provides  an 
explicit  analytical  expression  for  photon  distribution  function 
/(r,  s,  f),  as  a  function  of  position  r,  direction  of  light  s,  and 
time  t.  The  mean  position  and  the  half-width  at  half-maximum 
(HWHM)  height  of  the  distribution  are  always  exact.  In  this 
paper,  the  linear  forward  model  based  on  the  cumulant  solution 
is  described.  This  CFM  may  used  with  time-resolved,  contin¬ 
uous  wave  (CW),  and  frequency-domain  data,  which  are  much 
more  accurate  than  the  DA  models. 

To  obtain  a  3-D  image  one  needs  to  investigate  the  inverse 
algorithms.  For  clinical  applications,  this  requires  an  inversion 
technique,  that  is  computationally  fast,  and  stable  in  the 
presence  of  measurement  noise.  Recent  algorithms  to  solve  the 
inverse  problem  include  Newton’s  least-square-based  methods 
and  gradient-descent  methods  [l]-[5].  These  approaches  use 
an  iterative  procedure,  which  requires  a  long  computation 
time  to  solve  a  3-D  inverse  problem  with  large  unknowns 
(the  number  of  unknowns  is  the  number  of  voxels).  Further¬ 
more,  the  iterative  methods  can  not  ensure  that  the  result 
arrives  at  a  “global  minimum,”  and  does  not  converge  to  a 
“local  minimum,”  which  is  not  a  true  image  of  the  medium. 
The  application  of  Fourier  transform,  which  has  been  called 
“diffraction  tomography,”  can  greatly  reduce  computation 
time.  Matson  et  al.  [10]  and  Li  et  al.  [11]  have  developed  the 
diffraction  optical  tomographic  methods  to  realize  fast  image 
reconstruction.  However,  their  algorithms  are  limited  to  the  use 
of  a  single  light  source  with  a  two-dimensional  (2-D)  plane 
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of  detectors.  This  type  of  experimental  setup  acquires  only 
a  set  of  2-D  data  using  CW  or  frequency-modulated  light, 
that  is  not  enough  for  a  3-D  image  reconstruction.  Recently, 
Schotland  and  Markel  developed  inverse  inversion  algorithms 
using  diffusion  tomography  [12]-[14]  based  on  the  analytical 
form  of  the  Green’s  function  of  frequency-domain  diffusive 
waves,  and  point-like  absorbers  and  scatterers.  Data  obtained 
by  multiple  sources  with  multiple  detectors  in  parallel  slab 
geometry  are  used  in  these  approaches. 

A  fast  hybrid-dual-Fourier  (HDF)  algorithm,  which  uses  mul¬ 
tiple  sources  and  multiple  detectors  in  parallel  slab  geometry, 
is  described  in  this  paper  for  reconstruction  of  a  3-D  image  of 
an  inhomogeneous  medium.  This  approach  uses  a  general  2-D 
translation  invariance  of  the  Green’s  function  in  a  homogeneous 
background  slab  medium,  suitable  for  forward  models  based  on 
solution  of  RTE,  and  various  other  forward  models,  in  CW,  fre¬ 
quency-domain,  and  time-resolved  measurements.  This  inverse 
algorithm  runs  fast.  It  is  shown  that  a  3-D  image  of  a  turbid 
medium  (for  example,  divided  into  32  x  32  x  20  =  20480 
voxels)  can  be  reconstructed  in  1-2  min  using  a  personal  com¬ 
puter.  This  algorithm  can  produce  stable  images  in  presence  of 
relatively  strong  noises. 

The  forward  model  and  the  inverse  algorithm  discussed  in 
the  following  can  also  be  applied  for  image  reconstruction  in  a 
cloudy  environment  for  military  use. 

This  paper  is  organized  as  follows.  Section  II  presents  the 
analytical  solution  of  RTE,  based  on  a  cumulant  expansion, 
in  an  infinite  uniform  medium  and  shows  the  photon  distribu¬ 
tion  function  computed  using  the  cumulant  analytical  solution. 
Section  III  describes  the  forward  models  based  on  the  ana¬ 
lytical  solution  of  RTE,  considering  the  slab  geometry,  and  a 
weak  heterogeneity  using  a  perturbative  method.  Section  IV 
describes  the  HDF  inverse  algorithm  for  a  reconstruction  of 
a  3-D  image  of  an  inhomogeneous  medium.  The  3-D  image 
using  this  algorithm  is  shown.  A  discussion  is  presented  in 
Section  V. 

II.  Analytical  Cumulant  Solution  of  RTE 

The  photon  propagation  in  a  medium  is  described  by  the 
photon  distribution  function,  /(r,  s,  t),  as  a  function  of  time 
ty  position  r,  and  direction  s.  The  mathematical  equation 
governing  photon  propagation  is  the  well-known  radiative 
transfer  equation 

5J(r,  Syt)/dt-\-cs-  Vr/(r,  s,  t)  -f  /[Xa(r)/(r,  s,  t) 

=  j  -P(s,  s',  r)[/(r,  s',  t)  -  /(r,  s,i)]  ds' 

+  6{t  -  io)S{s  -  so)S{t  -  0)  (1) 

where  the  fundamental  parameters  are  the  scattering  rate 
/is(r)  =  cpcTs,  the  absorption  rate  /Xa(r)  =  cpcTa,  and  the 
differential  angular  scattering  rate  jUs(r)P(s,  s',  r),  where 
cTo  and  cTg  are  the  absorption  and  scattering  cross  sections 
respectively,  p  is  density  of  scatterers,  and  c  is  the  speed  of  light 
in  the  medium.  In  a  uniform  infinite  medium,  these  parameters 
are  position  independent. 


When  the  phase  function  depends  only  on  the  scattering 
angle,  we  can  expand  the  phase  function  in  Legendre  polyno¬ 
mials  with  constant  coefficients 

^(s,  s')  =  ^  ^  a(P,[cos(s  •  s')].  (2) 

Recently,  we  have  developed  a  new  approach  to  obtain  an 
analytical  solution  of  RTE,  based  on  a  cumulant  expansion,  in 
an  infinite  uniform  medium,  with  an  arbitrary  phase  function 
P{s,  s')  [8],  [9]. 

We  briefly  review  the  concept  of  “cumulant”  in  a  one-di¬ 
mensional  (1-D)  case.  Consider  a  random  variable  x,  with  a 
probability  distribution  function  /(x).  Instead  of  using  /(x) 
to  describe  the  distribution,  we  define  the  nth  moment  of 
X,  {x”)  =  /  x”/(.x)  ^^x,  and  correspondingly  the  nth  cumulant 
{x’")c  defined  by  cxp(^^^i(x")c(f0''/?^0  =  (exp(zfx))  = 

cumulant  (x}c  is  the  mean  posi¬ 
tion  of  X.  The  second  cumulant  {x^)c  represents  the  HWHM 
of  the  distribution.  The  higher  cumulants  are  related  to  the 
detailed  shape  of  the  distribution.  For  example,  (x^)c  describes 
the  skewness  or  asymmetry  of  the  distribution,  and  {x^)c 
describes  the  “kurtosis”  of  the  distribution,  that  is  the  extent  to 
which  it  differs  from  the  standard  bell  shape  associated  with  the 
normal  distribution  function.  The  cumulants,  hence,  describe 
the  distribution  in  an  intrinsic  way  by  subtracting  off  the  effects 
of  all  lower  order  moments.  In  3-D  case,  the  first  cumulant  has 
three  components,  the  second  cumulant  has  six  components, 
and  so  on. 

We  derived  an  explicit  algebraic  expression  of  spatial  cumu¬ 
lants  at  any  angle  and  any  time  that  is  exact  up  to  an  arbitrarily 
high  order  n  [9].  This  means  the  distribution  function  /(r,  s,  t) 
can  be  computed  to  any  desired  accuracy.  At  the  second  order, 
n  =  2,  an  analytic,  hence,  useful  explicit  expression  for  dis¬ 
tribution  function  /(r,  s,  t)  is  obtained  [8].  This  distribution  is 
Gaussian  in  position,  which  is  accurate  at  later  times,  but  only 
provides  the  exact  mean  position  and  the  exact  HWHM  at  early 
times.  A  weakness  of  the  second  order  cumulant  solution  is  that 
photons  at  the  front  edge  of  Gaussian  distribution  travel  faster 
than  light  speed,  thus  violate  causality,  though  to  a  much  less 
extent  than  that  in  the  DA. 

Fig.  1  compares  /(r,  s,  t)  obtained  from  the  analytical  cumu¬ 
lant  solution  and  the  Monte  Carlo  (MC)  simulation.  In  order  to 
reduce  the  statistical  deviation  to  an  acceptable  level,  10^  events 
are  counted  in  the  MC  simulation.  The  figure  shows  that  the 
solid  curve  (the  tenth-order  cumulant  solution)  is  located  in  the 
middle  of  data  obtained  by  the  MC  simulation.  The  solution  for 
CW  case  can  be  obtained  by  an  integration  of  /(r,  s,  t)  over 
time  t.  It  is  shown  that  even  second  order  cumulant  solution  (the 
dotted  curve)  can  provide  an  accurate  CW  solution,  because  this 
solution  ensures  that  the  mean  position  and  the  HWHM  of  dis¬ 
tribution  are  always  exact. 

The  plots  in  Fig.  1  indicate  that  a  strong  anisotropic  angular 
distribution  still  exists  at  z  ~  6  kr  {hr  is  the  transport  mean  free 
path)  from  the  source.  The  DA  is  only  valid  when  the  angular 
distribution  is  nearly  isotropic.  The  dominate  5  wave  distribu¬ 
tion  N{vy  ^)/47r  computed  using  the  diffusion  model  (the  thick 
dotted  curve)  has  a  large  discrepancy  with  the  MC  result. 
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Fig.  1.  Distribution  function  /(r,  s,  t)  in  an  infinite  uniform  scattering 
medium  as  a  function  of  time  t,  using  Henyey-Grccnstcin  phase  function 
with  g  —  0.9.  TTie  detector  is  located  at  i?  =  6  /tr  =  60/^  from  the  source 
front  along  direction  of  incident  light,  and  the  direction  is  along  the  incident 
direction.  The  solid  curve  is  computed  from  approximation  up  to  tenth  order  of 
cumulant;  the  dotted  curve  is  computed  from  approximation  up  to  the  second 
order  of  cumulant,  the  discrete  dots  are  from  the  MC  simulation;  the  curve  of 
thick  dots  is  from  the  DA,  N{x, 


where 

G  =  cQx^{-iXat)IF{s,  So,  t) 
fig)  =  [exp(f/i)  -  l]/g 
Ai  =  {l/4iT)exp{-git)  (6) 

Ty  is  obtained  by  replacing  cos  0  in  (5.2)  by  sin  0. 

The  HWHM  (second  cumulant)  is  expressed  as 

Bap{s,  t)  =  -  rlry2  (7) 


with 


(3)  ,  (^  +  1)^  p(4)' 
2/  +  3 


(8.1) 


^xx,yy  —  ^  ^  2  ■'4ir*i(C0S  0) 


Ki  - 1)  p.(i) 


{I  +  1){1  +  2)  (2)  l{l  -  1)  p,(3) 


2/ +  3 
{I  -f  !)(/  +  2)  (4) 

+  2/ +  3 . 


±Y,Ia,P^^\cosO) 


■  cos(20) 


^  pO)  I  ^  p(2) 


21-1 


The  second-order  analytical  cumulant  solution  is  given  by  [8] 
F(s,  So,  f)  1 


/(r,  s,t)  = 

where 
F(s,  So,  t) 


(47r)3/2  (detB)i/2 


•  exp 


-]{B-%0{r-r^Ur-n0 


(3) 


2/  +  1 

=  exp{-Hat)  ^  e\p{-git)Pi[cos{s  •  sq)] 

=  expi-Hat)  exp(-s/t)  ^yj„(s)y,;,(so). 


(4) 


In  (4),  g,  =  Ms[l  -  a,/{2l  +  1)],  YUO,  <P)  =  - 

(cos  6)  exp(2m0),  where  P/’^^(cos  6)  is 
the  associated  Legendre  function,  and  yim(s)  are  spherical  har¬ 
monics  normalized  to  An/ {21  H-  1). 

In  (3),  the  mean  position  of  the  distribution  (first  cumulant), 
when  the  source  is  located  at  tq  =  0  and  the  incident  direction 
is  along  z,  is  given  by 

r^(s,  l)=G^yl;P,(cosS) 

■[H+l)f{9i-0i+i)  +  lf{9i-gi-i)]  (5.1) 
r-'(s,  i)=G.^A,P/')(cos^) 

I 

■  cos 4>\f{gi  -gi-i)-  f{gi  -  (7/+1)]  (5.2') 


1 

21  -  1 


21-1-3^' 
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(8.2) 


where  (+)  corresponds  to  Axx  and  (—)  corresponds  to  Ayy. 
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f: 


(4) 


(8.4) 


Ayz  is  obtained  by  replacing  cos0  in  (8.4)  by  sin0.  In 
(8.1)-^8.4)  Fp“^^  are  given  by 

=  [/(<?/  -  91-2)  -  f{9i  -  9i-i)]/{9i-i  -  91-2)  (9.1) 

=  [f{9i  -  91^2)  -  f{9i  -  9i+i)]l{9w  -  91+2)  (9.2) 

E\^'^  =  \f{9i  -  m-i)  -  t]l{9i  -  9i-i)  (9.3) 

eI'^^  =  \f{9i  -  9i+i)  -  t]/i9i  -  9i+i)-  (9.4) 


Fig,  2(a)  and  (b)  shows  the  light  distribution  as  a  function  of 
time  at  different  receiving  angles  in  an  infinite  uniform  medium, 
computed  by  the  second  cumulant  solution,  where  detector  is 
located,  separately,  at  5  hr  [Fig-  2(a)]  and  15  hr  [Fig.  2(b)]  from 
the  source  in  the  incident  direction  of  the  source.  Fig.  2  shows 
the  existence  of  the  strong  anisotropy  of  the  light  distribution  at 
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t(ps) 

(a) 


t(ps) 


with  the  mean  position 

=  c[l  -exp{-git)]/gi,  (11) 


The  corresponding  time-dependent  diffusion  coefficients  are 


f  3^1  -  g2 


3i  Ui  9i{9i-92) 


[1  -  exp(-git)] 


[1  -  exp(-52^)] 


92(91  -  92) 

-^(1  -exp(-5ii)f  I 


29 


(12) 


D-s-j:  =  Dy 


+ 


92 


(13) 


As  shown  in  (1 1)-(13),  the  mean  position  of  the  distribution 
is  moving,  and  the  diffusion  coefficients  are  time  dependent.  At 
i  >  0,  the  mean  position  of  the  photon  density  moves  along 
z  direction  with  speed  c,  and  the  diffusion  coefficients  tend  to 
zero,  this  result  presents  a  clear  picture  of  near  ballistic  mo¬ 
tion.  As  time  increases,  the  mean  position  motion  slows  down, 
and  the  diffusion  coefficients  increase  from  zero.  This  stage  of 
photon  migration  is  often  called  a  snakelike  mode.  At  long  time, 
(10)  tends  to  the  center-moved  (1  Itr)  diffusion  model  with  the 
diffusion  coefficient  /tr/3. 


III.  Forward  Model  Based  on  the 
CUMULANT  Solution  of  RTE 


(b) 

Fig.  2.  Light  distribution  in  an  infinite  unifonn  medium  as  a  function  of  time 
at  different  received  angle,  using  second  cumulant  solution  of  radiative  transfer 
equation,  where  detector  is  located,  separately,  at  10  [Fig.  2(a)]  and  30  mm 
[Fig.  2(b)]  from  the  source  in  the  incident  direction.  The  parameters  for  this 
calculation  arc:  —  2  mm,  /«  =  300  mm,  the  phase  function  is  computed 

using  Mie  theory  for  polystyrene  spheres  with  diameter  cf  =  1.1 1  in  water 
and  the  wavelength  of  laser  source  A  =  625  nm,  which  gives  the  g-factor 
g  =  0-926. 


5  Ur  from  the  source  and  the  modest  anisotropy  at  a  distance  of 
15  Ur-  These  types  of  distributions  have  been  demonstrated  by 
time-resolved  experiments  [15], 

One  advantage  of  using  the  above  analytical  solution  of  RTE 
is  that  the  distribution  function  can  be  computed  very  fast.  The 
associated  Legendre  functions  can  be  accurately  computed 
using  recurrence  relations.  It  takes  only  a  minute  to  compute 
10^  data  of  /(r,  s,  t)  on  a  personal  computer. 

The  corresponding  solution  in  the  frequency-domain 
/(r,  s,  a;)  can  be  obtained  by  making  a  Fourier  transform 
/dt  exp(-za;f)/(r,  s,  i).  The  CW  solution  is  obtained  by 
taking  u;  =  0. 

The  photon  density  N{r,  t)  of  the  second  cumulant  solution 
is  given  by 


iV(r,  t)  = 


exp 


{z  - 


•exp 


{x^  -f 


4Da:xCt 


AD..ct 


exp(-^a0  (10) 


The  linear  forward  models  for  scattering  media  are  built  in 
following  three  steps:  1)  computation  of  a  background  Green’s 
function  in  an  infinite  uniform  medium;  2)  extension  of  this 
Green’s  function  to  slab  geometry;  and  3)  computation  of  the 
weight  function  using  a  perturbative  method.  These  steps  have 
been  applied  in  building  the  linear  forward  models  under  DA 
[2].  We  use  these  steps  as  well,  but  our  approach  is  based  on  the 
cumulant  solution  of  RTE,  rather  than  the  solution  of  the  diffu¬ 
sion  equation. 

We  use  the  second-order  cumulant  solution  for  computing  a 
background  Green’s  function  in  an  infinite  uniform  medium, 
since  it  is  easy  to  use  the  explicit  expressions  in  (3)-<9),  that 
avoid  complicated  computations  of  higher  order  cumulants.  The 
second  order  cumulant  solution  is  accurate  at  later  times,  but 
only  provides  the  correct  mean  position  and  the  correct  HWHM 
at  early  times.  We  notice  that  the  width  of  the  distribution  at 
early  times  could  be  smaller  than  the  size  of  a  voxel,  the  average 
over  the  distributions  at  different  points  in  a  voxel  smears  the  de¬ 
tail  shape  of  the  distribution.  In  the  CW  or  frequency-  domain 
cases,  the  shape  of  the  distribution  is  further  smeared  by  integra¬ 
tion  over  time  t.  Therefore,  the  second-order  cumulant  solution 
can  be  a  reasonable  approximation  in  building  forward  models 
based  on  the  RTE. 

Since  a  detector  usually  collects  emergent  light  within  a  wide 
range  of  angle  of  different  directions,  it  is  convenient  to  compute 
the  Green’s  function  related  to  a  detector  using  photon  density 
t)  (10)~(13),  where  Vd  is  the  position  of  detector. 
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Fig.  3.  Schematic  diagram  shows  how  to  extend  the  cumulant  solution  of  RTE 
from  an  infinite  medium  to  a  semi-infinite  medium. 


It  is  essential  to  include  the  boundary  effect  in  the  solution  of 
the  RTE  when  photons  are  injected  into  and  spread  out  from  a 
finite  sized  medium.  A  proper  extension  of  the  cumulant  solu¬ 
tion  to  slab  geometry  is  an  essential  step  for  building  a  forward 
model. 

A  boundary  condition  is  applied  based  on  the  following  phys¬ 
ical  consideration.  At  early  times,  the  center  of  photon  distribu¬ 
tion  injected  into  medium,  moves  forward  into  medium.  Then, 
the  distribution  spreads  out  from  the  moving  center  with  dif¬ 
fusion  coefficients  that  gradually  increase  from  zero.  At  early 
times,  the  number  of  photons  leaking  out  of  the  boundary  is 
negligible  compared  to  the  total  number  of  the  incident  photons. 
The  boundary  condition  plays  a  role  at  later  times,  when  there 
are  many  photons  leaking  out  of  the  boundary. 

The  approach  known  as  an  approximate  “extrapolated” 
boundary  condition  [16],  extrapolates  the  boundary  by  a 
distance  ^  =  otltn  the  extrapolation  length,  beyond  the  real 
boundaries  with  a  ~  0.7,  at  which  the  photon  density  vanishes. 

To  apply  this  boundary  condition  for  the  cumulant  solution  in 
a  semi-infinite  geometry  a,  a  virtual  negative  source  Sy  is  added 
to  the  original  source  5,  as  shown  in  Fig.  3.  During  the  early 
period,  the  solution  of  the  RTE  in  an  infinite  uniform  medium 
automatically  satisfied  the  boundary  condition  because  the  den¬ 
sity  is  near  zero  at  the  boundary,  and  the  virtual  source  does  not 
play  a  role.  After  a  time  of  approximately  4  Itr/c,  the  center  of 
photon  density,  C,  has  moved  and  stopped  at  a  position  1  /tr 
from  the  original  source  S  and  the  center  from  virtual  source, 
Cv,  has  moved  in  a  similar  way.  Then,  the  arrangement  shown 
in  Fig.  3,  produces  a  cancellation  of  contributions  to  the  photon 
density  from  the  original  source  and  the  virtual  source  on  the 
extrapolated  boundary. 

Fig.  4  shows  that  the  time-resolved  backscattered  photon  dis¬ 
tribution  in  a  semi-infinite  medium  on  the  =  0  surface,  with 
the  source-detector  distance  1  Ur,  obtained  using  the  second- 
order  cumulant  approximation  and  the  extrapolated  boundary 


Fig.  4.  Backscattcrcd  photon  distribution  /(r,  s  =  — z,  t)  emerging  from 
plane  surface  of  a  semi-infinite  turbid  medium,  as  a  function  of  time,  with  the 
source-detector  distance  1  /tr  on  the  surface  z  ^  0  plane.  The  pulse  source 
is  located  at  ^  =  0,  incident  along  z  direction.  The  extrapolated  boundary 
condition  is  used.  The  solid  curve  is  obtained  from  cumulant  approximation 
(CA),  up  to  the  second  cumulant.  The  dashed  curve  is  from  DA.  The  cross  points 
are  obtained  from  MC  simulation. 

condition,  which  agrees  with  the  MC  simulation  much  better 
than  that  of  the  DA. 

For  extending  to  the  slab  geometry,  adding  a  series  of  pairs  of 
virtual  “image”  sources  at  both  sides  of  slab  is  a  good  approxi¬ 
mation  for  satisfaction  of  the  extrapolated  boundary  conditions 
on  both  sides  of  a  slab  [17]. 

The  heterogeneous  structure  of  a  highly  scattering  turbid 
medium  can  be  characterized  by  the  following  optical  param¬ 
eters:  the  scattering  rate  /is(r),  the  absorption  rate  //a(r),  and 
the  differential  angular  scattering  rate  fis{r)P{s,  s\  r). 

A  perturbation  method  is  used  which  takes  the  photon  dis¬ 
tribution  function  in  a  uniform  background  slab  medium  as  the 
zero-order  approximation.  The  change  of  the  photon  distribu¬ 
tion  function  originates  from  the  change  of  optical  parameters 
compared  to  that  in  the  uniform  background  slab  medium.  The 
change  of  scattering  and  absorption  parameters  are  defined  as 
follows: 

Afisir)  =Ms{r) 

A/ia(r)  =Ma(r)  - 

A[msP](s,  s',  r)  =/Xs(r)P(s,  s',  r),  s')  (14) 

where  the  quantities  with  super  index  (0)  are  the  optical  pa¬ 
rameters  in  a  uniform  background  slab  medium.  By  expanding 
A[/is^]  (s,  s',  r)  in  Legendre  polynomials,  we  obtain 

A[/isP](s,  s',  r) 

=  ^  ;^[A/i,(r)a(“>  +  Mi“)Aa,(r)]P,lcos(s  •  s')]  (15) 

with  Aao(r)  =  0,  since  ao  always  equals  1.  The  physical 
meaning  is  that  the  scattering  parameters  have  no  effect  on  the 
s  (I  =  0)  component. 
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Making  a  perturbation  expansion  of  (1)  to  the  first-order  Bom 
approximation,  the  change  in  the  photon  distribution  is  given  by 

A/(rd,  Sd,  firs,  Ss) 

dr  j  c/s'/^°^(rd,  ss,t-t'\r,  s') 

•  jy  A[/isP](s,  s',  r)7W(r,  s,<'  |r„  s,}ds 

-  [AM,(r)  +  A^la{T)]I^°\T,  s',  t'  I  rs,  s,)|  (16) 

where  A7(rd,  Sd,  i  |  Ss)  is  the  change  in  the  light  intensity 
received  by  a  detector  located  at  rd,  along  the  direction  Sd,  and 
at  time  which  is  injected  from  a  source  located  at  r^,  along 
a  direction  of  s.,,  at  time  f  =  0.  “Change”  refers  to  the  dif¬ 
ference  in  intensity  compared  to  that  received  by  the  same  de¬ 
tector,  from  the  same  source,  when  light  passes  through  a  uni¬ 
form  background  slab  medium.  The  term /^^^(r2,  S2,  f  |ri,  si) 
is  the  intensity  of  light,  calculated  using  the  cumulant  solution 
of  RTE,  at  r2  along  the  direction  S2  and  at  time  t,  when  light  is 
injected  from  a  position  ri  along  a  direction  of  si  at  time  i  =  0 
migrating  in  a  uniform  background  slab  medium. 

The  background  Green’s  functions  in  (16),  obtained  by  cu¬ 
mulant  solution,  are  expanded  in  spherical  harmonics 

/(®^(r,  s,t'  |r„  Ss)  =  r^,  s^,  t')Ytmis), 

'  I,  m 


7(°)(rd,  sa,  t-t'\T,s)  =  sa,  t  -  <')y,:„(s). 

(17) 


The  spherical  transform  is  performed  using  a  fast  Fourier  trans¬ 
form  for  the  integral  over  0,  and  a  Clenshaw-Curtis  quadrature 
for  the  integral  over  0, 

Using  the  orthogonality  relation  of  the  spherical  function  and 
the  addition  theorem:  ^/m(s)V/*,(s')  =  Pz[cos(s  ♦  s')],  the 

analytical  integration  over  s  and  s'  in  (1 6)  can  be  performed.  For 
time  resolved  data,  the  contribution  from  an  absorbing  object 
located  at  is  given  by 


A7(rd,  Sd,  Is,  Ss,  t  |rfc) 

""AMa(rfc)6yfc^  ‘^^'2  (2m) 

■'^Alm(Tk,  Is,  Ss,  t')Ci^{Tk,  Td,  Sd,  t  -  t')  (18) 


where  6Vk  is  the  volume  of  kth  voxel,  and  L  is  the  cutoff  value 
in  the  Legendre  expansion  in  (18).  The  contribution  from  a  scat¬ 
tering  object  located  at  Vk  is  given  by 


A/(rcj,  Ts,  Ss,  t  j  r^) 


rt  L 

I 

1  =  1 


47r 


A/is(rfc)  1  ~ 


{21  -f  1) 

21  +  1 


..(0)  Aai(rfc) 
2/  +  1 


’  ^  ^  Alrnij^kj  ^s,  Ss,  t  (ffc,  I'd,  Sd ,  ^  ^  ) .  (1  9) 


For  frequency  domain  (or  CW)  data,  the  contribution  from  an 
absorbing  object  located  at  rjt  is  given  by 

A7(rd,  Sd,  Is,  Ss,  wlr/t) 


L 


= -Atia{rk)SVkY^ 

/=0 


47r 

2/  +  1 


,„(rfc,  fs,  Ss,  u))c*^{Tk,  rd,  Sd,  w)  (20) 

rrt 


and  the  contribution  from  a  scattering  object  located  at  rfc  is 
given  by 


A7(rd,  Sd,  Is,  Ss,  u;|rt) 


/=1 


47r 

2/TT 


AAt.,(rk)  1 


..(0)  Aaz(rfc) 


2/  -h  1  y  "  2/  -M 

•  Ai,n{Tk,  Fs,  Ss,  a;)4„(rfc,  rd,  Sd,  w). 


(21) 


Comparing  (1 8)-(2 1)  with  the  corresponding  weight  function 
commonly  used  in  the  DA,  [1],  [2]  only  s  wave  (/  =  0)  for  ab¬ 
sorptive  objects,  and  only  p  wave  (/  =  1)  for  scattering  objects 
are  considered  in  the  diffusion  forward  models.  Besides,  even 
for  s  wave  and  p  wave,  the  diffusive  solution  is  incorrect  when 
voxels  are  located  near  the  source,  as  discussed  before. 

The  previous  formulae  allow  simulating  the  background 
Green’s  function  and  the  change  of  optical  parameters  in 
detail.  They  are  also  applicable  to  the  cases  where  only  a 
few  parameters  of  the  medium  are  known,  similar  to  that 
for  the  diffusion  forward  model.  When  only  \  and 

p-factor  for  an  uniform  background  medium  are  given,  the 
Henyey-Greenstein  phase  function  [18]  is  widely  adopted  as 
an  approximate  phase  function 

1  1  - 

i^(cos^)  =  —  ^|^^2E2pcosl9)^ 

=  ^^(2Z  +  l)<7'P,(cos^).  (22) 


Although  (22)  uses  a  single  parameter  p-factor  to  describe  a 
phase  function,  this  description  is  much  better  than  that  used 
in  the  DA,  which  implies  a  phase  function  linear  in  cos  0, 

If  Aai  (r)  in  (2 1 ),  which  represent  the  change  of  the  phase 
function,  is  not  considered,  two  optical  parameters  being  im¬ 
aged  are  Apa(r)  and  Aps(r).  The  reduced  scattering  coeffi¬ 
cient  A/is(l  -  is  directly  related  to  AD  (change  of  the 

diffusion  coefficient)  used  in  the  DA  models.  The  CFM,  hence, 
can  be  applied  to  the  experimental  data  in  a  similar  fashion  as 
that  for  the  DA  models,  to  obtain  images  of  the  optical  parame¬ 
ters.  In  the  CFM,  however,  all  contributions  from  higher  spher¬ 
ical  waves  are  properly  included. 

The  most  time  consuming  part  in  computation  of  CFM  using 
the  previous  formulae  is  to  build  a  database  of  Aim  and  CL- 
Once  it  is  built  for  a  uniform  background  medium,  the  database 
can  be  applied  for  imaging  of  various  heterogeneity  cases.  In 
parallel  geometry,  Aim  is  a  function  of  (xk  —  Xs,  Vk  “  2/s)  due 
to  the  2-D  translation  invariance.  Since  position  of  source  2:^  and 
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incident  direction  Ss  are  fixed,  only  a  3-D  (xjk-Xs,  yk-ys,  ^k) 
database  is  required.  When  Ss  is  taken  along  2;  direction  (light 
is  injected  perpendicular  to  surface),  the  scale  of  database  is  re¬ 
duced  to  2-D  due  to  the  z  axis  symmetry.  Photons  from  different 
directions  in  a  wide  solid  angle  are  received  by  a  detector,  as  dis¬ 
cussed  before,  photon  density  N{vd  —  r,  t)  is  used  for  com¬ 
puting  the  Green’s  function  associated  with  detectors,  which  is 
independent  of  and  can  be  computed  much  easily.  The 
database  can  be  built  in  a  reasonable  computation  time  because 
the  distribution  function  (r2,  S2,  t\ri,  Si)  can  be  rapidly 
calculated  using  the  analytical  expressions. 


IV.  Fast  3-D  HDF  Inverse  Algorithm 

We  now  outline  an  inverse  algorithm  to  quickly  reconstruct 
image  of  a  medium  from  acquired  measurements  using  the 
above  CFM.  The  model,  neglecting  the  irrelevant  parameters, 
can  be  briefly  written  as 


T'lfd,  fs,  Zs) 


j  dYdzW{Yd  -  f,  Fs  -  f,  z,  Zd,  z) 


(23) 


where  R  =  (f,  z)  is  the  position  of  a  voxel  inside  turbid 
medium;  f  is  (x,  y)  coordinates;  Rg  —  (fs,  Zg)  is  the  position 
of  a  source;  and  Rd  ==  (fd,  zj)  is  the  position  of  a  detector. 
In  (23),  y(rd,  fg,  z^,  Zg)  is  the  measured  change  in  light 
intensity  received  by  a  detector  at  Rd  from  a  point  source  at  Rg. 
A’(r,  z)  is  the  change  of  the  optical  parameters  inside  turbid 
medium.  The  weight  function  IV(rd  —  r ,  fg  —  r,  z,  z^,  Zg) 
is  a  function  of  Fd  ~  r  and  fg  —  r  on  (x,  y)  plane,  because  of 
parallel  geometry,  assuming  an  infinite  sized  area,  and  the  2-D 
translation  invariance  of  the  Green’s  function  in  a  background 
homogeneous  slab.  Here,  the  special  form  of  the  weight 
function  is  not  relevant;  the  weight  function  can  be  calculated 
by  the  CFM  or  the  DA  models,  using  with  CW,  frequency, 
or  time-resolved  data.  This  approach  is  general  and  can  also 
be  used  for  inverse  problems  of  nonoptical  measurements  in 
parallel  geometries. 

A  light  source  scans  through  a  2-D  array.  Transmitted  or 
backscattered  light  signals  emerging  from  the  medium  are 
detected  using  a  2-D  array  of  detectors,  such  as  a  charge-cou¬ 
pled  device  (CCD)  camera  (or  time-gated  CCD  camera  in 
the  time  resolved  case).  Each  illumination  of  the  light  source 
provides  a  set  of  2-D  data  on  the  2-D  detector  array.  For  CW 
or  frequency-modulated  light  source,  this  arrangement  can 
produce  a  set  of  2-D  x  2-D  =  4-D  data  in  a  relatively  short 
acquisition  time,  because  a  CCD  camera  produces  2-D  data 
of  the  detectors  at  different  positions  simultaneously.  When 
time-resolved  or  modulation  at  multiple  frequencies  are  ap¬ 
plied,  a  set  of  five-dimensional  (5-D)  data  can  be  acquired.  The 
inverse  problems  of  3-D  imaging,  hence,  are  over-determined, 
which  is  necessary  for  obtaining  an  accurate  3-D  image. 

When  the  translation  invariance  is  satisfied,  the  Fourier 
transform  approach  is  a  powerful  technique  to  achieve  a  fast 
inversion.  In  the  Fourier  space,  the  convolution  of  W  and  X 
becomes  a  product  of  W  and  X,  and  the  weight  matrix  W 
becomes  diagonal.  Hence,  inversion  can  be  performed  much 
faster.  Using  this  concept  in  the  case  of  multiple  sources  and 


multiple  detectors  in  parallel  geometries  a  dual  2-D  Fourier 
transform  /  dr^  is  performed  on  (23),  to  obtain 


j  f^,  z^,  Zd) 

~  I  J  I  £f(rd  —  f 

.  W{t,  ~  r,  rd  -  r,  z„  za,  z)e'(^^-"^^>^X(f,  z) 


which  leads  to 


y(qd,  qs,  Zd,  Zs)  =  /  dzW{qd,  Qs,  2:,  Zd,  2^5)^ (qd  H-  qs,  z) 

.  .  .  (24) 

where  Y,  X,  and  W  are  change  in  light  intensity,  change  in 
optical  parameters,  and  the  weight  function  in  the  Fourier  space, 
respectively. 

A  similar  form  of  this  dual  Fourier  transform  has  been  derived 
by  Markel  and  Schotland  [13],  [14]  in  a  frequency-domain  dif¬ 
fusion  model. 

Equation  (24)  seems  difficult  to  be  used  for  performing  the 
inverse  reconstruction  because  of  the  argument  mismatch  (qj  + 
Qs)  in  X  and  (qs,  qu)  in  y  and  W.  This  difficulty  occurs  be¬ 
cause  the  weight  function  in  (23)  is  related  to  three  positions: 
Fd,  Fs,  and  F.  To  remove  this  complexity,  the  following  linear 
hybrid  transform  is  introduced; 

u  =  qd  +  qs 

v=qd-qs.  (25) 

This  results  in  the  HDF  formula 

y(u,  V,  Zd,  Zs)  -  J  dzW{n,  V,  z,  Zd,  Zs)X{n,  z)  (26) 

where  Y,X,  and  W  are,  respectively,  Y,  X,  and  W  as  functions 
of  u  and  v. 

While  (25)  is  a  relatively  simple  expression,  it  is  essential  to 
properly  realize  this  hybrid  transform  in  discrete  lattices  of  the 
Fourier  space.  A  procedure  to  quickly  perform  this  transform 
from  (qd,  qs)  coordinates  to  new  (u,  v)  coordinates,  separately, 
for  X  and  y  components,  is  explained  in  Fig.  5  using  an  example 
of  a  6  X  6  lattice.  The  maximum  value  of  u  is  taken  as  the 
maximum  value  of  qd  or  qs,  not  the  maximum  value  of  qd  +  qs- 
The  periodic  property  of  lattice  in  the  Fourier  space  is  used,  for 
example,  y(u  =  2,  v  =  4)  -  y(qs  =  3,  qa  =  5),  This 
procedure  builds  a  one-to-one  correspondence  between  lattices 
in  the  two  coordinate  systems.  Fig.  5  shows  that  Y  and  W  at 
each  node  [circle  in  Fig.  5]  in  (u,  v)  coordinates  are  directly 
mapped  from  Y  and  Wy  respectively,  at  the  corresponding  node 
in  (qd,  qa)  coordinates  without  any  algebraic  manipulation. 

In  (26),  a  common  2-D  Fourier  argument  u  appears 
in  y,  X,  and  W,  For  each  value  of  u,  (26)  leads  to  an 
over-determined  1-D  problem  for  inverse  reconstruction: 

=  /  dzW{v,  z)X{z).  In  order  to  perform  fast  in¬ 
version,  we  invert  the  normal  form  of  the  forward  model: 
YW'^  =  [W'^W]X  for  each  u,  where  [W'^W]  is  a  M  x  M 
matrix,  with  M  the  number  of  layers  in  z  direction.  The  original 
W  in  (23)  is  a  matrix  with  a  large  dimension.  The  inverse 
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Fig.  5.  Example  of  a  6  X  6  lattice  for  explaining  the  linear  hybrid  transform 
from  (<ict,  fL)  coordinates  to  (ti,  v)  coordinates. 


problem  now  is  simplified  to  invert  many  (number  of  discrete 
value  of  u)  matrices,  each  with  a  small  dimension  M.  The  latter 
problem  is  much  more  computationally  efficient  compared 
to  the  original  problem  of  (23).  Once  X(u,  z)  are  obtained 
for  all  u,  a  2-D  inverse  Fourier  transform  produces  X(r,  z), 
which  is  the  3-D  image  of  optical  parameters  of  the  medium. 
Markel  and  Schotland  use  different  procedures  for  inversion.  In 
[13],  a  Fourier-Laplace  inversion  is  applied,  hence,  an  analytic 
continuation  of  measured  data  to  the  complex  plane  is  required 
for  the  inverse  Laplace  transform.  In  [14],  an  inverse  procedure 
is  performed  in  an  argument  space,  similar  to  variables  v  here. 
Since  v  include  2-D  variables,  inversion  in  v  space  could  take 
longer  time  than  that  of  inversion  in  2:  space. 

As  discussed  previously,  matrices  W  and  [I'F^W^]  for  each 
u  can  be  calculated  in  advance  for  a  uniform  background  slab 
medium.  Assuming  that  a  group  of  experimental  data  has  been 
acquired,  the  following  steps  are  taken  to  produce  a  3-D  image 
of  the  medium: 

1)  obtain  “change”  of  intensities,  Y  (fa,  fs,  z^j  Zg),  by  sub¬ 
tracting  the  intensity  for  a  uniform  background  medium 
from  the  measured  intensity; 

2)  extend  the  (x,  y)  area  and  padding  zeros,  to  overcome  the 
wraparound  problem  in  discrete  convolutions  [19]; 

3)  perform  a  dual  2-D  fast  Fourier  transform  (FFT)  of 

rs»  ^d,  ^a)  in  the  extended  area  to  produce 
y(qdj  qs,  ^a); 

4)  determine  y(u,  v,  2^,  Zg)  for  each  u,  using  a  mapping 
procedure  explained  in  Fig^  5; 

5)  invert  YW'^  =  \W'^W]X  for  each  u,  which  is  an  in¬ 
verse  problem  involving  a  M  x  M  matrix,  with  M  the 
number  of  layers  along  2:  direction.  Proper  regularization 
according  to  noise  level  needs  to  be  taken  into  account. 
Regularization  will  be  discussed  later  in  the  paper; 

6)  perform  an  inverse  2-D  FFT  on  X(u,  z)  to  produce 

X(f,  2). 


Fig.  6.  3-D  image  reconstructed  using  hybrid  dual  Fourier  tomography.  Two 
absorbing  objects,  each  with  the  volume  3x3x2  mm^,  arc  located  inside  a 
turbid  medium  with  volume  96  X  96  x  40  mm^  divided  into  32  x  32  X  20 
voxels.  The  first  one  is  located  at  position  labeled  (10,  10,  10)  with  absorption 
difference  A/i^  =  0.01  mm“L  The  second  one  is  located  at  position  labeled 
(20,  20,  15)  with  absorption  difference  A/io  =  0.007  mm~L  A  CW  light 
source  incident  perpendicular  to  the  z,  =  0  plane  is  scanned  through  a  2D  32  x 
32  array  at  the  plane,  with  each  pixel  3  mm  x  3  mm.  A  same  sized  2D  array  of 
detectors  is  located  at  Zd  plane  (transmission  geometry).  The  simulated  data  arc 
produced  with  noise  5%.  A  linear  scale  of  color  bar  from  the  maximum  value  to 
minimum  value  of  A^i^  is  used.  The  numbers  labels  the  z  layers  counting  form 
source  to  the  detector,  layers  arc  separated  by  2  mm. 

Our  computational  experiments  show  it  takes  only  1-2  min 
on  a  personal  computer  to  perform  an  inverse  reconstruction  of 
a  3-D  image  of  a  medium  with  a  large  number  of  voxels  (for 
example,  32  x  32  x  20  voxels)  using  this  HDF  algorithm. 

To  demonstrate  our  concept  of  HDF  tomography  in  3-D 
image  reconstruction,  an  example  using  simulated  CW  data  is 
presented.  A  slab  turbid  medium,  with  a  transport  mean-free 
path  Itr  =  1  mm,  absorption  length  la  —  300  mm,  and 
thickness  Zd  =  40  mm,  is  divided  into  20  layers.  A  CW  light 
source,  injected  perpendicular  to  the  Zg  =  0  plane,  scans  by 
a  2-D  32  X  32  array  on  the  plane,  with  each  pixel  3x3  mm. 
A  2-D  array  of  detectors  with  the  same  spacing  is  located  at 
Zd  plane  (transmission  geometry).  The  medium,  is  divided  into 
32  X  32  X  20  voxels,  each  of  dimension  3x3x2  mm^. 
Two  absorbing  objects  are  located  in  the  medium,  each  with  a 
volume  3x3x2  mm^.  The  first  one  located  at  (10,  10,  10) 
has  an  absorption  difference  of  A/Xa  =  0.01  mm~^  with  the 
background.  The  second  one  is  located  at  (20,  20,  15)  with  an 
absorption  difference  of  0.007  mm“^.  The  simulated  data  with 
noise  level  of  5%  are  obtained  using  the  CFM.  The  tomographic 
images  are  shown  in  Fig.  6.  As  shown,  the  central  positions  of 
3-D  image  of  the  objects  are  correct,  located  at  a  voxel  (10, 10, 
10)  with  dark  color,  and  a  voxel  (20,  20,  15)  with  gray  level. 
The  resolution  of  image  is  about  ~6  mm  in  the  transverse 
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(x,  y)  plane  and  ~10  mm  along  2:  direction.  In  general,  the 
axial  resolution  (along  2:  direction)  is  poorer  than  the  lateral 
resolution  [on  the  (t,  y)  plane].  In  transmission  geometry,  two 
Green’s  functions  in  the  weight  function  compensate  each  other 
when  the  2:  position  of  the  object  changes,  that  leads  to  a  poor 
sensitivity  of  the  measured  photon  intensity  to  the  z  position 
of  the  object.  The  shapes  of  3-D  image  of  two  objects  are 
ellipsoids  with  longer  axis  along  the  2:  direction.  The  absorption 
difference  has  the  maximum  value  at  the  center  of  ellipsoid, 
and  decays  gradually  with  increase  distance  from  the  center. 

A  cutoff  in  discrete  lattices  of  qs  and  naturally  introduces 
a  kind  of  regularization.  This  regularization  is  very  effective. 
Initial  tests  show  that  even  adding  30%  of  fluctuations  on  sim¬ 
ulated  data  of  y(fd,  Fs,  2^,  2.,,),  an  image  similar  to  that  shown 
in  Fig.  6  is  still  reconstructed.  The  reason  for  this  is  that  noises 
come  from  fluctuations  at  different  source  and  detector  posi¬ 
tions,  which  are  mainly  the  high-frequency  components  of  qs 
and  qa.  A  cutoff  in  qs  and  q^  naturally  eliminates  these  high 
frequency  noises,  such  that  a  stable  image,  especially  in  (x,  y) 
plane,  can  be  reconstructed  in  a  strong  noise  level. 

However,  the  inverse  problem  is  still  ill-posed,  because 
contribution  to  the  change  of  intensity  from  a  small  voxel 
deeply  inside  medium  is  weak,  and  is  not  sensitive  to  its  2 
position  in  transmission  case.  A  regularization  procedure  on 
inversion  of  71^^  -  [W'^WjX  is  still  needed.  The  standard 
Tikhonov  regularization  approach  [20]  is  applied  and  L-curve 
[21],  [22]  method  is  used  for  determining  the  best  regularization 
parameters. 

This  fast  inverse  algorithm  produces  a  3-D  image  in  a  linear 
image  regime.  For  nonlinear  image  reconstruction  procedure, 
the  reconstructed  3-D  image  provides  a  good  initial  profile  for 
further  refining  the  3-D  image  taking  the  nonlinear  effects  into 
consideration. 

The  HDF  inversion  method  can  be  extended  to  a  cylindrical 
geometry,  with  an  arbitrary  shape  of  the  (x,  y)  cross  section, 
for  3-D  image  reconstruction.  In  this  geometry,  an  algorithm 
using  a  single  Fourier  inversion  has  been  developed  [23].  This 
algorithm  is  limited  to  the  case  that  the  sources  and  the  detectors 
are  located  on  a  plane  with  same  2  coordinates.  The  hybrid-dual- 
Fourier  inverse  approach  in  cylindrical  geometry  removes  this 
restriction,  so  more  data  can  be  acquired  for  3-D  tomography. 
The  linear  forward  model  in  cylindrical  geometry  is  given  by 

Fg,  2d,  Zs} 

drdzWivd,  Fg,  F;2d  -  2,  2^  -  z)X{t,  2)  (27) 

where  VF(Fd,  Fg,  F;2d  -  2,  2^  -  2)  is  the  weight  function,  a 
function  of  Zd  -  2  and  Zg  —  2  due  to  the  1-D  translation  invari¬ 
ance  of  the  Green’s  function  in  a  homogeneous  background 
medium  in  cylindrical  geometry  (assuming  infinite  2  length). 
We  make  a  dual  1-D  (along  2  direction)  Fourier  transform 
J dzddzse^<^^'^e^^^‘  on  (27)  to  obtain 

^(qd,  qs,  Fd,  Fg)  =  /  dzW{q^,  qs,  F,  Fd,  Fs);^(qd  +  qg,  F) 

.  .  .  (28) 
where  7,  X,  and  W  are  the  Fourier  space  quantities  corre¬ 
sponding  that  in  (27). 


The  (1-D)  linear  hybrid  coordinate  transforms,  u  =  qd  +  qg, 
and  v  =:  qd  -  qs,  for  (28)  leads  to 

y(u,  V,  Tj,  r,)  =:  J d?W(u,  V,  fd,  r,;r)X(u,  r)  (29) 

where  7,  X,  and  W  are,  respectively,  7,  X,  and  W  as 
functions  of  u  and  v.  For  each  value  of  u,  (29)  is  an  over 
determined  2-D  problem  for  inverse  reconstruction,  namely,  to 
determine  a  2-D  unknown  value  of  X(u,  r)  from  known  3-D 
data  of  7(u,  v,  Fd,  Fg)  for  each  u.  This  3-D-2-D  determination 
enhances  the  accuracy  of  3-D  image  compared  to  2-D-2-D 
determination  in  the  single-Fourier  transform  inversion.  After 
X{u,  r)  are  obtained  for  all  u,  a  1-D  inverse  Fourier  transform 
produces  the  image  X{r,  2). 

V.  Discussion 

As  shown  in  (19)  and  (21),  there  is  no  contribution  from  5 
wave  to  the  weight  function  for  a  scattering  object.  This  result 
reflects  a  fact  that  no  scattering  effect  exists  for  an  isotropic  an¬ 
gular  distribution.  In  the  regions  far  from  sources,  the  weight 
function  contributed  from  scattering  objects  is  small  because 
there  is  no  contribution  from  the  dominant  s  wave,  as  shown  in 
many  results  based  on  the  diffusion  models  [1H5].  This  non¬ 
sensitivity  of  signals  to  the  scattering  objects  deep  inside  the 
medium  should  be  considered  in  optical  tomography.  A  pure 
isotropic  distribution  is  never  achieved,  otherwise,  there  will  be 
no  flux  in  any  directions.  In  the  diffusive  model,  a  small  p  wave, 
-  (3/47r)Ds  -VN,  exists  which  maintains  the  photons  diffusing 
to  the  regions  with  fewer  photons.  The  factor  —  VW  represents 
this  effect.  However,  this  expression  is  valid  only  in  the  regions 
where  the  p  wave  is  much  smaller  than  s  wave,  (l/47r)iV,  and 
does  not  correctly  describe  the  early  photon  propagation  near 
sources.  Since  only  the  weight  function  for  scattering  objects 
close  to  sources  plays  an  important  role,  but  it  was  estimated 
using  the  formula  valid  in  regions  far  from  sources,  substantial 
error  introduced  in  the  diffusion  forward  model  for  scattering 
objects  is  crucial,  considering  >  pa  in  tissue. 

For  the  weight  function  of  absorbing  objects,  contributions 
from  all  spherical  components,  including  s  wave,  are  given  in 
(18)  and  (20).  In  commonly  used  diffusion  formula,  the  con¬ 
tribution  from  p  wave  was  neglected.  The  diffusion  coefficient 
originally  derived  in  the  DA  is  D  =  1/(3^^  4-  p^),  that  leads  to 
AD  =  The  contribution  from  p  wave  to 

the  weight  function  for  absorbing  objects,  hence,  should  exist. 
However,  in  the  later  diffusion  models,  AD  is  assigned  only 
for  scattering  objects  and  only  5  wave  for  absorbing  objects  is 
taken.  Equations  ( 1 8)  and  (20)  provide  a  quantitative  estimation 
of  weight  function  for  absorbing  objects  in  regions  close  to  the 
source,  as  well  as  far  from  the  source. 

The  CFM  and  the  HDF  inverse  algorithm  need  further 
improvements  in  the  following  aspects.  Further  improvement 
should  be  considered  without  significantly  increasing  com¬ 
plexity  in  computation.  First,  the  second  cumulant  solution  is 
not  accurate  in  the  detailed  shape  of  the  distribution,  especially, 
the  front  edge  in  the  Gaussian  distribution  violates  causality. 
An  empirical  distribution,  which  keeps  the  exact  value  of  the 
first  and  second  cumulants,  while  satisfies  the  causality,  can  be 
designed  to  replace  the  Gaussian  distribution. 
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Second,  the  boundary  condition  is  approximate.  When  a 
more  accurate  distribution  /(r,  s,  t)  at  early  time  is  needed, 
the  boundary  condition  for  a  semi-infinite  geometry  should  be 

I{x,  y,  2  =  0;  e,  (^,  t)  =  0,  if  cos  >  0.  (30) 

This  type  of  the  boundary  condition  was  studied  by  Domke  [24] 
for  the  steady  state  case.  The  solution  is  represented  as  a  super¬ 
position  of  a  solution  describing  a  transport  problem  in  an  infi¬ 
nite  medium,  and  a  Fredholm  integral  term,  which  corrects  this 
solution  for  the  appropriate  half-space  boundary  condition.  This 
approach  may  be  used  for  further  development  of  the  boundary 
problem. 

Third,  to  consider  the  nonlinear  effects,  in  (16)  should 
be  replaced  by  the  Green’s  function  in  a  real  heterogeneous 
medium.  Among  the  high-order  perturbative  corrections  of 
the  Green’s  function,  the  ‘‘self-energy”  diagram,  which  counts 
photon  round  trips  through  a  position  up  to  infinite  times, 
plays  an  important  role.  Gandjbakhche  e(  al  [25]  studied  this 
effect  using  a  random  walk  model.  We  find  that  a  renormal¬ 
ization  procedure  for  this  nonlinear  effect  can  be  performed 
after  image  is  obtained  using  a  linear  inversion  process.  This 
renormalization  procedure  can  recover  the  optimal  value  of  the 
optical  parameters  and  can  improve  the  resolution  of  image. 
The  detailed  results  of  the  renormalization  will  be  published 
elsewhere. 

The  translation  invariance  is  valid  for  the  parallel  geometry 
assuming  that  the  (x,  y)  area  is  infinite.  We  suppose  that  this 
assumption  of  the  infinite  area  is  reasonable.  How  much  error 
arises  due  to  the  finite  area  of  a  sample  will  be  studied  in  details. 

Use  of  the  simulated  data  mainly  tests  the  validity  of  the 
inverse  algorithm,  does  not  test  accuracy  of  the  forward  model. 
Experimental  data  from  phantoms  and  in  vivo  measurements 
in  human  body  will  be  performed  for  further  testing  of  our 
approach. 
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Abstract 

The  powerlaw  patterns  in  Mie  scattering  (the  normalized  light  intensity  I(0)/I{0)  vs.  the  dimensionless  qR  where 
q  =  4nr'  sin  |  is  the  magnitude  of  the  wave  vector  transfer  at  the  scattering  angle  0  for  wavelength  A,  and  R  is  the  radius 
of  the  nonabsorbing  sphere  with  a  relative  refractive  index  m  >  1)  are  analyzed  using  the  geometrical  optics  approx¬ 
imation  for  particles  of  a  large  size  parameter.  The  powerlaw  regime  is  shown  to  be  present  only  in  Mie  scat¬ 

tering  of  soft  particles.  The  {qR)~^  powerlaw  regime  occurs  at  the  scattering  angles  of  the  /?  =  1  geometrical  ray 
(refracted  without  internal  reflections)  from  the  portion  of  the  incident  beam  with  an  incidence  angle  around  7t/4  upon 
the  particle.  The  {qR)'^  powerlaw  regimes  from  particles  sharing  one  common  relative  refractive  index  but  differing  in 
size  parameters  are  collinear.  Simple  analytical  expressions  are  derived  to  describe  these  powerlaw  regimes  of  Mie 
scattering. 

©  2003  Elsevier  B.V.  All  rights  reserved. 

PACS:  3.80.-fr;  42.25.Fx;  78.35+c 
Keywords:  Light  scattering;  Mie  scattering 


The  study  of  light  scattering  by  small  particles  is 
important  in  noninvasive  characterization  of  small 
particles  and  radiative  transfer  in  turbid  media 
including  atmosphere,  marine  environment  and 
tissues  (see,  for  example,  van  de  Hulst’s  classic 
work[l]  and  the  recent  review  volume  edited  by 
Mishchenko,  Hovenier  and  Travis  [2]).  Light 
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scattering  from  a  sphere  of  arbitrary  size  and  re¬ 
fractive  index  (Mie  scattering),  one  of  a  few  exactly 
solvable  cases,  was  derived  in  1908  [3].  This  exact 
solution  was  given  in  the  form  of  a  slow  con¬ 
verging  partial-wave  series  involving  complex 
functions.  The  physical  meaning  and  interpreta¬ 
tion  of  the  Mie  scattering  was  itself  of  lasting  in¬ 
terests  [4-7], 

Recently,  work  on  powerlaw  regimes  in 
Mie  scattering  was  obtained  by  Sorensen  and 
Fischbach  [8]  when  plotting  the  normalized  light 
intensity  7(0) //(O)  versus  the  dimensionless 
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parameter  gR,  where  ^  =  is  the  mag¬ 

nitude  of  the  wave  vector  transfer  at  the  scatter¬ 
ing  angle  6  for  wavelength  A  after  ignoring  the 
interference  ripple  structure.  These  patterns  were 
attributed  to  the  structure  factor  of  the  illumi¬ 
nated  portion  of  the  scattering  object.  The  Fou¬ 
rier  transform  of  the  illuminated  annular  shell  for 
a  sphere  of  radius  R  was  used  to  explain  the 
emerging  [gR)^,  powerlaw  re¬ 

gimes.  This  approach  ignores  the  extra  phase  shift 
incurred  to  the  light  when  it  passes  through  the 
particle.  The  implicit  assumption  made  there  [8] 
that  the  phase  shift  due  to  the  nonunity  refractive 
index  of  the  particle  is  negligible  is  valid  for 
scattering  of  X-rays  [9],  but  it  is  inappropriate  for 
optical  light  scattering  from  particles  with  a  large 
phase  shift. 

In  this  paper,  we  analyze  patterns  in  Mie 
scattering  using  the  geometrical  optics  approxi¬ 
mation  (GOA)  for  particles  of  a  large  size  pa¬ 
rameter  and  derive  simple  analytical  expressions 
to  describe  these  powerlaw  regimes  of  Mie  scat¬ 
tering. 

The  increase  of  the  phase  shift  accompanies 
with  the  increase  of  the  size  parameter.  The  GOA 
becomes  a  viable  one  when  the  size  parameter 
X  =  /:/?>  1  where  k  ~  2n?r^  is  the  wave  number. 
The  scattering  amplitude  of  light  is  composed  of 
a  diffraction  light  component  and  reflected  and 
refracted  rays  when  the  contribution  from  surface 
waves  can  be  neglected  [10,1].  The  diffraction 
peak  is  highly  concentrated  around  the  exact 
forward  direction  within  angle  A9  ^  l/x.  By  ig¬ 
noring  the  interference  ripple  structure,  the  Mie 
scattering  in  near  forward  scattering  directions 
and  outside  the  diffraction  peak  for  soft  particles 
whose  refractive  index  m  is  close  to  unity 
(|/M  -  1|  <  1)  and  for  dense  particles  {m  ^  1.5)  in 
GOA  is  dominated  by  the  p  ~  1  geometrical  ray 
(refracted  without  internal  reflections)  [10-12]. 
The  {gR)"^  and  {gR)~^  powerlaw  regimes  can  be 
recovered  from  the  asymptotic  behavior  of  the 
contribution  from  this  p—\  geometrical  ray  [see 
Fig.  1]. 

The  magnitude  of  the  scattering  amplitude  of 
perpendicular  polarization  with  respect  to  the 
scattering  plane  which  was  examined  in  [8]  (the 
parallel  polarization  case  can  be  analyzed  in  a 


Fig.  1.  Geometrical  rays  scattered  by  a  sphere. 


similar  fashion)  is  dominated  by  the  contribution 
Ar^  from  the  p  =  1  geometrical  ray  [10] 

-  1 )'( 1+  ^  + 1^) 

(1) 


and  the  scattering  angle  9  is  expressed  as 


.  9  f(Vm2  +  (w2-  l)?-l) 

Sm  -  =  — - 7- - - 

2  m(H-/2) 


in  terms  of  the  incidence  angle  t  where  t  =  tan  t. 

For  soft  particles,  Eqs.  (1)  and  (2)  can  be  ex¬ 
panded  in  powers  of  /z  =  w  -  1 


=  ^t  +  0{l?), 


(3) 

(4) 


yielding 


X  p  X  p 

2  pi  A-  sin^  I  2  sin^  | 


(5) 


under  the  condition  /i  sin  -  ^  y/Ji,  The  first 
approximate  form  in  (5)  appeared  in  the  classic 
book  of  van  de  Hulst  [1,  p.  222]. 

When  sin  the  incidence  angle  r  is  around 
/  =  tanT=l,  Eqs.  (1)  and  (2)  can  also  be  ex¬ 
panded  about  this  point  t  ==  1 
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Since  the  refractive  index  m  ~  1, 


(6) 


(V) 


_ w^  +  1 _ , 

+1) 

for  soft  particles,  and  the  two  values  inside  the 
brackets  of  (6)  and  (7)  are  close  even  for  dense 
particles,  we  find 


^  \/2x  {s/lrn^  —  1)^^^ 

'"■  “sin|  (v/2^nrf+  1)5/2- 


(8) 


Fig.  2.  The  ratio  (Ak^  sin|)/(/t/ji  sinf),^,  versus  /  -  tanr.  The 
curves  obtained  for  particles  of  different  refractive  indices  in¬ 
tersect  with  the  90%  line  at  /  0.7  and  ?  ~  1.6  where  the  value 

of  the  incidence  angle  is  t  =  0.27i  and  t  =  0.3n,  respectively. 


The  range  of  the  incidence  angle  t  over  which  the 
above  expression  (8)  is  valid  can  be  examined  by 
plotting  {An,  sin|)/(y4«,  sinf),^,  versus  /  =  tanr  [see 
Fig.  2].  "The  range  of  the  incidence  angle  is  approx¬ 
imately  0.2n  <  T  <  0.371  for  soft  and  dense  particles 
(1  <m  <  2)  when  a  10%  relative  error  is  allowed. 
The  corresponding  scattering  angle  range  is  given  by 

-  1  -  \/2  .  d  V35m^  -  25  -  v/IO 

3m  ^*'"2^  Im 

(9) 

Thus,  we  can  write  the  normalized  scattered  light 
intensity  as 


(b)  qR 

Fig.  3.  The  powerlaw  regimes  of  a  sphere  of  refractive  index; 
(a)  m  =  1.01;  (b)  w  =  1.50.  The  size  parameter  is  x  =  384. 
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from  Eqs.  (5)  and  (8)  in  GOA  where  we  have  used 
the  fact  S\  (0)  =  x^/2  and  qR  =  2x  sin  f.  This  shows 
that  the  {qRy^  powerlaw  regime  exists  in  Mie 
scattering  of  both  soft  and  dense  particles  while 
the  {qRy^  powerlaw  regime  only  appears  in  soft 


particles.  The  {qR)~^  powerlaw  regimes  are  col- 
linear  while  the  {qRy  powerlaw  regime  are  not 
for  particles  sharing  one  common  relative  refrac¬ 
tive  index  but  differing  in  size  parameters.  For 
example,  the  {qRy^  powerlaw  regime  is  found 
within  5  <  qR  <  12,  followed  by  the  {qRy"^  pow¬ 
erlaw  regime  within  8  <  <  80,  for  a  soft 

sphere  of  a  refractive  index  1.01;  only  the 
{qRy^  powerlaw  regime  is  observed  within 
165  <  qR  <  300  for  the  dense  sphere  of  a  refractive 
index  m  =  1.50  (see  Fig.  3).  The  size  parameter  of 
the  sphere  in  both  cases  is  x  =  384. 

Fig.  4  demonstrated  the  powerlaw  regimes  for  a 
nonabsorbing  sphere  with  refractive  indices 


Fig.  4.  The  normalized  Mie  scattering  curves  plotted  versus  qR  for  spheres  of  refractive  index:  (a)  w  =  1.01;  (b)  w  =  1 .05;  (c)  m  =  1.5. 
The  {qR)~^  and  {qR)~^  powerlaw  regimes  of  Mie  scattering  given  by  Eq.  (10)  are  also  plotted. 
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m  =  1.01,  m  =  1.05  and  w  =  1.5,  respectively.  The 
[qRy^  powerlaw  regime  exists  only  in  soft  parti¬ 
cles  and  disappears  in  dense  particles.  Its  trend 
agrees  well  with  our  simple  expression  (10)  [see  the 
long-dash,  dash,  dot  and  dash-dot  lines  for  parti¬ 
cles  of  increasing  size  parameters  in  Fig.  4(a)  and 
(b)].  This  agreement  is  better  for  larger  and  softer 
particles.  On  the  other  hand,  the  {qRY^  powerlaw 
regime  exists  in  both  soft  and  dense  particles.  This 
regime  occurs  at  a  larger  value  of  qR  for  particles 
of  a  larger  size  parameter  and  is  broader  for  denser 
particles.  The  {qR)~^  powerlaw  regimes  of  particles 
of  a  common  refractive  index  but  of  different  size 
parameters  coincide  on  one  straight  line  (the  solid 
lines  in  Fig.  4(a)-(c)). 

In  conclusion,  we  have  analyzed  the  powerlaw 
patterns  in  Mie  scattering  using  the  geometrical 
optics  approximation  for  particles  of  a  large  size 
parameter.  The  [qRY^^  powerlaw  regime  is  shown  to 
be  present  only  in  soft  particles.  The  {qR)"^  pow¬ 
erlaw  regime  occurs  at  the  scattering  angles  of  the 
p  —  1  geometrical  ray  (refracted  without  internal 
reflections)  from  the  portion  of  the  incident  beam 
with  an  incidence  angle  around  7t/4  (from  0.2tc  to 
0.3tc  within  a  10%  relative  error  of  the  scattering 
amplitude)  upon  the  particle.  The  {qRY^  powerlaw 
regimes  from  particles  sharing  one  common  relative 
refractive  index  but  differing  in  size  parameters  are 
collinear.  The  {qRY^  and  {qRY^  powerlaw  regimes 
of  Mie  scattering  are  well  captured  by  the  simple 
analytical  expressions  given  in  Eq.  (10). 
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ABSTRACT 

We  report  on  the  effect  of  the  nonlinear  multiple  passage  on  optical  imaging  of  an  absorption  inhomo¬ 
geneity  of  finite  size  deep  inside  a  turbid  medium  based  on  a  curaulant  solution  to  radiative  transfer. 
An  analytical  expression  for  the  nonlinear  correction  factor  is  derived.  Comp^son  to  Monte  Carlo 
simulations  reveals  an  excellent  agreement.  The  implication  on  optical  imaging  is  discussed. 

Keywords;  nonlinear  correction,  multiple  passage,  radiative  transfer,  optical  imaging 

1.  INTRODUCTION 

The  principle  of  optical  imaging  of  turbid  media  (such  as  tissues)  is  to  locate  and  reconstruct  the  optical 
properties  (absorption  and  scattering  coefficients)  of  embedded  inhoraogeneities  (such  as  tumor)  in  the 
hope  of  identification  by  inverting  the  difference  in  time-resolved  or  frequency-modulated  photon  trans¬ 
mittance  due  to  the  presence  of  the  inhomogeneities  through  either  iterative  or  noniterative  methods. 
The  key  quantity  involved  is  the  weight  function  which  quantifies  the  influence  on  the  detected  signal 
due  to  the  change  of  the  optical  parameters  of  the  medium.  The  diffusion  approximation  to  radiative 
transfer  provides  an  adequate  model  for  the  weight  function  (or  Jacobian)  for  a  small  and  weak  ab¬ 
sorption  inhomogeneity  far  away  from  both  the  source  and  the  detector.  However,  the  weight  function 
predicted  by  the  linear  perturbation  approaches  is  no  longer  valid  when  the  absorption  strength  is  not 
small.  ^  This  can  be  attributed  to  the  multiple  passage  of  a  photon  through  one  single  abnormal  site. 

The  change  of  the  light  intensity  A/  at  the  detector  due  to  the  presence  of  an  absorption  site  at 
r  from  a  modulated  point  source  at  Tg  is  expressed  as 

AJ  =  -<i/iarG(rd,wlr)G(r,a;!r,)  (1) 

to  the  first  order  of  Bom  approximation  where  J/Xo  is  the  excess  absorption  of  the  absorption  site  whose 
volume  is  V,  cv  is  the  modulation  frequency  of  light,  and  G  is  the  propagator  of  photon  migration  in 
the  background  medium.  Here,  the  Green’s  function  G(r2,a;|ri),  in  general,  depends  on  the  detail  of 
light  scattering  inside  the  medium,  and  the  incident  and  outgoing  directions  of  light. 

When  the  absorption  strength  is  not  small  {SfiaV  ^  1),  photon  loss  due  to  multiple  passage  of  the 
absorption  site  is  appreciable  8ind  can  not  be  ignored.  The  expression  for  A/  in  Eq.  (1)  needs  to  be 
modified  to  include  the  contributions  from  multiple  visits  of  the  site  by  the  photon.  Fig.  (1)  illusti'ates 
the  most  important  correction  (a  “self-energy”  correction)  which  takes  into  account  the  repeated  visits 
made  by  a  photon  to  the  site  up  to  an  infinite  times. 
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Figure  1.  Self-energy  correction  to  the  multiple  passage  effect,  on  light  absorption. 


Assuming  that  the  center  of  the  absorption  site  is  located  at  f  and  far  away  from  both  the  source 
and  the  detector,  the  change  of  the  detected  light,  A/,  is  now  given  by 


n5=0 

-G(rrf,  wir)  ^ 


(2) 


where 

((j;  =  ^  fyf^  <'(*'2,  tij|ri)d®r2d®ri  (3) 

is  the  self-propagator  which  describes  the  probability  that  a  photon  revisits  the  volume  V  of  size  R. 
Here  (?(rd,u;|f)  and  G(f,a/|ra)  are  well  modelled  by  the  center-moved  diffusion  model  as  long  as  the 
separations  jr^  —  f  | ,  |ra  —  r|  »  i{  where  It  is  the  transport  mean  free  path  of  light  in  the  medium. ^ 
However,  the  diffusion  Green’s  function  can  not  be  used  in  Eq.  (3)  to  evaluate  N^i{{uj',R)  where  ri  is 
in  the  proximity  of  ra.  By  comparing  Eq.  (2)  to  Eq,  (1),  the  nonlinear  multiple  passgage  effect  of  an 
absorption  site  can  be  summarized  by  the  nonlinear  correction  factor  [1  +  Naeifiu)',  ii)V’5/iQ(f)]"\  This 
factor  serves  as  a  universal  measure  of  the  nonlinear  multiple  passage  effect  as  long  as  the  absorption 
site  is  far  away  from  both  the  source  and  the  detector  and  its  size  is  much  smaller  than  its  distance  to 
both  the  source  and  the  detector. 


In  this  article,  we  will  derive  an  analytical  expression  for  the  self-propagator  to  understand  the 
nonlinear  multiple  passage  effect  on  light  absorption  using  our  cumulant  solution  to  radiative  transfer. 
The  nonlinear  correction  factor  [l  +  (w|  i?)V5jUa(r)]“^  of  our  result  is  shown  to  be  in  an  exceUcnt 

agreement  with  the  Monte  Carlo  simulations  for  continuous  wave  light. 


2.  THEORY 

To  take  into  account  the  higher  order  contributions  from  the  absorption  inhomogeneity,  the  behavior 
of  the  photon  migration  in  a  short  distance  must  be  considered.  Although  the  photon  distribution  is 
almost  isotropic  at  an  absorption  site  deep  inside  the  medium,  the  diffusion  approximation  is  still  not 
appropriate  here.  The  separation  between  the  two  points  ri  and  ra  within  the  volume  in  Eq.  (3)  is 
smaU.  The  photon  propagator  which  repr^ents  the  probability  that  a  photon  propagates 
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from  position  n  with  propagation  direction  s  to  position  r2  in  time  t,  wheti  rz  is  in  the  proximity  of 
Ti,  is  governed  by  the  radiative  transfer  equation  rather  than  the  difiusion  equation. 

Recently  we  have  shown  that  the  propagation  of  photon  inside  a  turbid  medium  (the  radiative 
transfer  equation)  can  be  solved  analytically  using  a  cmnulant  expansion  of  the  photon  distribution 
function.®  The  propagation  of  photon  was  found  to  transform  from  an  initial  ballistic  motion  at  early 
time  and  then  gradually  to  a  center-adjusted  diffusion  at  later  time.  The  propagator  of  photon  density 
(the  Green’s  function)  in  an  infinite  uniform  medium  is  ^ven  by’* 


'I''”  — 4D(t)t '  - 

ignoring  the  small  difference  in  the  diffusion  coefficient  along  different  directions  where  the  absorption 
coefficient  is  Ha,  the  time-dependent  diffusion  coefficient  is 

D{t)  =  ^  I  ^  -  [1  -  exp(-ci/lt)]  -  I  [1  -  exp(-ct/4)]®|  (5) 


(r  -  rp  -  8oA(f))® 
AD{i)t 


A(t)  =  /*(!  -  exp(-ct//t)]  (6) 

is  the  average  center  of  photons  which  moves  with  speed  c  initially  and  approaches  the  transport  mean 
free  path  It  in  the  long  time  limit.  The  Green’s  function  for  parallel  geometries  can  be  obtained  by  the 
method  of  image  sources.® 

2.1.  Propagator  of  an  isotropic  point  source 

Let’s  now  consider  the  propagator  iV(r,  tjrcSo)  at  the  inhomogeneity  site  rp  =  0  (the  origin  of  space) 
deep  inside  the  medium.  The  photon  distribution  at  rp  is  almost  isotropic  but  is  anisotropic  scattering. 
The  effective  propagator  can  then  be  obtained  by  averaging  (4)  over  the  propagation  direction  sp  of 
light  over  the  47r  solid  angle,  and  is  given  by  [see  Appendix  A] 


^  I  ,Ps,N{r.t\ro.^)  ^ 


TA(t) 


:|exp  - 


(r-A«))- 


-  exp  “ 


(r  +  A(t)) 

AD(t)t 


This  reduces  to 


^cfr(»‘,  t) 


exp(-Mat) 


|«qj  1^- 


(r  -  hf 

iDoat 


-  (47r)3/a(I?^t)t/Mt  . . 

in  early  and  late  time  limits  where  Doo  =  itc/3. 

The  temporal  Fourier  transforms  of  the  asymptotic  equations  (8)  and  (9)  are  given  by 


{r  +  kf 

ADex^t 


for  t »  1 
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aad 

^eff(r,w)  =  g^j^^[exp(-K|r-Zt|)-exp(-«(r  +  it))]  (11) 

respectively,  where  k  b.  \/3(fia  —  iuj)fltc  whose  sign  is  chosen  with  a  nonnegative  real  part.  In  the  limit 
of  small  K  1,  Eq.  (11)  simplifies  to 


Urn  ATeffCr.u;) 


=1 


1 

47rDoor 


r  <  It 
r>lt 


(12) 


This  is  the  case,  for  example,  that  a  continuous  wave  propagates  in  a  nonabsorbing  medium.  The 
erromous  divergence  at  the  zero  separation  in  the  diffuse  Green’s  function 


G(r,w) 


exp(— Kr) 
AttDooT 


(13) 


is  removed  in  our  formulation  of  the  propagation  of  an  isotropic  point  source. 

The  asjmptotic  equation  (11)  from  the  late  time  limit  provides  a  good  approximation  for  Neff(r,u/) 
when  r>  It  [see  Fig.  (2)].  The  contribution  to  Neff(r,ui)  when  r  <  is  from  either  ballistic  or  diffusive 
photons,  hence  an  improvement  to  Eq.  (10)  can  be  made 


+  r<i,  (14) 

to  include  the  contribution  from  diffusive  photons.  The  effective  propagator  in  temporal  Fourier  space 
Areff(r,w  -  0)  and  its  asymptotic  behaviors  (10),  (11)  and  (14)  are  shown  in  Fig.  (2).  The  diffusion 
Green’s  function  has  a  huge  error  for  small  r. 


2.2.  Self  propagator  for  a  finite  volume 

For  an  absorption  site  of  a  finite  volume  V  deep  inside  the  medium,  say  a  sphere  of  radius  R  ^  L 
where  L  is  the  dimension  of  the  medium,  the  selfrpropagator  f^eelf(*;  R)  for  this  volume  which  denotes 
a  photon  revisits  the  site  in  time  t  is  written  as: 


where 


7o(r) 


Iviv  ”  ri| ,  t}dhid^T2 

1 

V  Jo 


(15) 


(16) 


is  the  characteristic  function  for  a  uniform  sphere.®"®  This  characteristic  function  has  a  form  of 

70(0  =  1  -  {S/4V)r  +  . . .  (17) 

for  an  arbitrary  particle  where  S  is  the  surface  area  of  the  particle.  This  self  propagator  (15)  for  a  finite 
volume  is  quite  different  from  the  self-propagator  of  a  point,  obtained  by  setting  r  =  0  in  (4)  or  (7), 


expi-fXgt)  Ajtf ' 

(47rI?(t)t)3/2  [  AD{t)t\  ’ 


(t  >  0). 


(18) 
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Figure  2.  The  effective  propagator  in  temporal  Fourier  space  =  0)  for  photon  migration  in  a  nonab- 

Eorbing  medium.  Its  approximations  by  (14)  when  r  <  It  and  by  Eq.  (11)  when  r  >  It  are  also  plotted.  Thn 
diffusion  Green’s  function  has  a  huge  error  for  small  r. 


See  Fig.  (3).  This  difference  comes  from  the  fact  that  Eq.  (15)  includes  the  contribution,  from  the 
ballistic  motion  of  the  photon  when  the  photon  Sies  across  the  site  while  Eq.  (18)  does  not  contain  this 
effect.  The  former  manifests  itself  in  Fig.  (3a)  as  the  linear  decay  of  hi  the  form  of  7o(ci) 

near  the  origin. 

The  self-propagator  in  temporal  Fourier  space  is  thus  obtained  by  a  temporal  Fourier  transform  of 
(15); 


iV'aelf  J?)  = 


A+OO  _ 

/  (t;  R)  exp  (iuit)  dt 

J0+ 

=  f  dr^Q{r)4:-KT^  /  ci«Nejr(r,  t)exp(fwt) 
V  Jo  Jo 

1 

~  Y  iVeff(r,u;)7o(r)47rr^dr. 


(19) 


The  lower  limit  of  integration  is  0^,  emphasizing  that  t  =  0  should  be  excluded  from  integration.  Note 
limf_*o+  t)  =  0  for  our  cumulant  photon  density  function.  This  is  not  the  case  for  the  diffusion 

Green's  function.  A  numerical  quadrature  is  generally  required  to  compute  this  self  propagator  (19).  A 
crude  estimation  of  ■^)  can  be  obtained  from  the  as5rmptotic  behavior  (11)  and  (14)  of  NestirfU)), 

i.e., 


1  rminl 

^yl  j;^exp 

+i  r  ' 

Vjo 


8iTDfx,rKtt 


{iu>  —  A*o)-J  'YQ{r)4Trr^dr  (20) 

[exp(— « |r  —  itl)  —  exp(— K(r  +  it))]  70  (r)47r7-''^  dr. 
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Figure  3.  The  self-propagators  for  a  finite  volume  and  a  point;  (a)  Nseifit’,  R)  and  (b)  iVcff(0,  t). 


This  reduces  to 


'SR  , 

TB  +  m 

384ft5+1601?il3-6CM?fl=+31? 

-SWR^hc 


R  <  lt/2 
R  >  hl2 


for  a  continuous  wave  propagating  inside  a  nonabsorbing  medium  (w  =  //„  =  «  =  0).  This  estimation 
turns  out  to  be  amazingly  good.  Fig.  (4)  plots  iV^jf  (u;  =  0;  R)  from  numerical  quadrature  and  the  crude 
estimation  (21). 


3.  RESULTS  AND  DISCUSSION 

The  multiple  passage  effect  due  to  the  absorption  site  can  now  be  computed  using  the  self-propagator 
Eq.  (19)  derived  here.  For  large  sites,  the  self-propagator  iVgeif  (w  =  0;  B)  increases  inverse  proportional 
to  its  size  (iVaeif  oc  R~^)  from  Eq.  (21);  hence  the  nonlinear  correction  factor  has  a  form  of 


dependent  on  the  area  of  the  absorption  site  for  large  R. 

Monte  Carlo  methods  have  been  extensively  used  in  simulation  of  photon  migration.®- We  perform 
Monte  Carlo  simulations  on  a  uniform  nonabsorbing  and  isotropic  scattering  slab  (the  anisotropic  factor 
of  scattering  g  *=  0).  The  units  of  length  of  time  are  chosen  such  that  the  mean  scattering  length 
la  -  l/fis  =  1  and  the  speed  of  light  c  =  1.  The  transport  mean  free  path  is  hence  It  =  1  and  the 
thickness  of  the  slab  is  assumed  L  =  80ft  ■  An  absorption  spherical  site  is  located  at  the  center  of  the 
slab  (0,0,Z,/2)  with  radius  R  whose  absorption  and  scattering  coefficients  are  ^a,2  =  S/Xa  =  0.01  and 
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Figure  4.  The  self  propagator  JV’seif(u;  =  0;fl)  and  its  estimator.  The  diffusion  self-propagator  for  continuous 
waves  is  also  plotted. 

Ma,2  =  tis  respectively.  The  photon  is  incident  at  the  origin  on  the  left  boundary  of  the  slab  z  =  0  in 
the  normal  direction  of  the  surface.  Each  photon  is  traced  until  it  escapes  the  slab  through  either  the 
left  or  the  right  boundary.  The  correlated  sampling  is  used  in  simulation  to  reduce  variance.  A  single 
simulation  is  used  to  compute  the  emitted  photon  deusily  Iq  for  the  uniform  background  (nonabsorption 
slab)  and  I  for  the  slab  with  the  absorption  site  present. 

The  nonlinear  correction  factor  [l -I- in  Eq.  (22)  can  be  extracted  from  the 
change  of  the  detected  light  intensity  due  to  the  presence  of  the  absorption  site  in  Monte  Carlo  simula¬ 
tions  according  to  Eq.  (2).  Fig.  (5)  plots  the  theoretical  nonlinear  correction  factor  and  that  from  Monte 
Carlo  simulations.  “Back”  and  “Forward”  denote  the  cases  where  light  emits  from  the  left  (z=0)  and 
the  right  {z  =  L)  boundaries,  respectively.  The  agreement  between  our  theoretical  result  and  Monte 
Carlo  simulations  is  excellent  except  for  extremely  small  sizes  of  inhomogeneities. 

Figs.  (6)  and  (7)  plot  the  nonlinear  correction  factor  versus  the  variation  of  the  modulation  frequency 
of  light  for  a  fixed  absorption  strength  and  versus  the  size  of  the  absorption  site  with  a  fixed  modulation 
frequency  of  light  respectively.  With  the  increase  of  the  modulation  frequency  of  light,  the  nonlinear 
correction  becomes  less  accentuated.  The  dependence  on  the  size  of  the  inhomogeneity  is  no  longer 
monotonic  for  modulated  light  while  the  nonlinear  correction  factor  decreases  monotonically  with  the 
increase  of  the  size  for  continuous  wave  light.  The  phase  delay  is  in  the  order  of  a  few  degrees  in  the 
cases  investigated. 

The  typical  value  of  the  absorption  coefficient  of  human  tissues  is  around  O.OOlps"^  while  the 
scattering  coefficient  is  about  Ips"^.  Hence  the  absorption  and  scattering  ratio  is  in  the  order  of  0.001. 
This  should  be  compared  to  our  results  listed  here  where  the  corresponding  ratio  is  0.01  and  one  order 
of  magnitude  stronger.  The  nonlinear  correction  factor  for  absorption  inhomogeneities  such  as  tumors 
in  human  tissues  is  not  appreciable  unless  the  size  of  the  inhomogeneity  is  R  bit  or  larger. 
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Figure  5.  The  nonlinear  correction  factor  from  the  theoretical  self-propagator  Eq.  (19)  and  Monte  Carlo  simu¬ 
lations.  “Back”  and  “Forward”  denote  light  emitting  from  the  left  (z=0)and  the  right  (z  =  t)  boundaries.  The 
excess  absorption  is  ==  0.01. 
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Figure  6.  The  nonlinear  correction  factor  versus  the  variation  of  the  modulation  frequency  of  light  The  size  of 
the  absorption  sphere  is  R  =  3lt.  The  excess  absorption  is  6/ia  =  0-01. 
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Figure  7.  The  nonlinear  correction  factor  versus  the  variation  of  the  size  of  the  absorption  site.  The  modulation 
frequency  of  light  is  w  =  0.1.  The  excess  absorption  is  6tia  =  0.01. 

In  conclusion,  we  have  derived  an  analytical  expression  for  the  nonlinear  correction  factor  which 
agrees  well  with  Monte  Carlo  simulations.  The  ^ect  of  the  nonlinear  multiple  passage  of  an  absorption 
site  on  optical  imaging  only  becomes  appreciable  when  the  size  of  the  inhomogeneity  is  5/t  or  larger  for 
human  tissues. 


APPENDIX  A.  DERIVATION  OF  NEFwiR,  T) 

The  spatial  Fourier  transform  of  (4)  is  given  by 

iV(k,f|ro,so)  =  j d®rexp(— 2k-r)fV(r,t|ro,so)  =exp(-fc^D(t)t-/iat-  tk-soA(t)^  .  (23) 

Hence,  the  effective  propagator  in  spatial  Fourier  space  at  ro  is  expressed  as 


Ar^(k,«)  d2solV(k,f|ro,so)  -  exp  -  iij) 


(24) 


by  averaging  (23)  over  the  propagation  direction  Sq  of  light  over  the  47r  solid  angle.  The  effective 
propagator  in  real  space  is  then  obtained  by  an  inverse  spatial  Fourier  transform  of  (24): 


^efT(r,  t) 


f  d^k 
J  (2vr)3 


exp(ik  •  r)  exp  ^—k^D{t)t  —  fiat^ 


sin(fcA(t)) 

A;A(t) 


exp(-/xat) 

(47r)3/2(D(t)i)V2rA(i) 

2exp(— //of) 
(47r)3/3(D(t)f)V2rA(<) 


(r  -  A(t)) 

r^  +  A{tf' 
W{t)t 


2) 


-  exp 


(r  +  m) 

4Z)(t)t 


sinh 


rA(t) 


(25) 
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characterization  of  absorptive  inhomogeneities 
embedded  in  turbid  media 
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Abstract:  Independent  component  analysis  of  the  scattered  wave  is  proposed  as  the  preproces¬ 
sor  to  characterization  of  absorptive  inhomogeneities  embedded  in  turbid  media.  Reconstruction 
results  with  simulated  and  experimental  data  will  be  presented  for  multiple  embedded  objects. 
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OCIS  codes:  (170.3010)  Image  reconstruction  techniques;  (290.7050)  Turbid  media;  (100.3190)  Inverse  problems 

Optical  biomedical  imaging  has  attracted  significant  interests  over  the  past  decade  because  of  its  potential  to 
noninvasively  probe  the  interior  of  turbid  media  such  as  human  breasts. [1,  2,  3]  The  foreign  objects  such  as 
tumors  inside  the  turbid  medium  distort  the  light  flux  surrounding  the  medium.  The  distortion  of  the  light 
flux  (the  scattered  wave)  can  be  used  to  reconstruct  the  location  and  the  optical  properties  of  the  foreign 
objects  via  modelling  of  the  light  propagation  inside  the  medium. [4]  In  the  linearized  scheme  of  inversion, 
the  perturbation  on  the  light  flux  due  to  the  inhomogeneity  is  determined  by  two  Green’s  functions:  the  first 
one  G(r,  Ts)  of  light  propagation  from  the  source  at  to  the  site  r  of  the  inhomogeneity,  and  the  second 
one  G'(rd,r)  of  light  propagation  from  the  site  of  inhomogeneity  to  the  detector  at  r^.  The  inverse  problem 
is  highly  ill-posed  because  of  the  smoothing  nature  of  the  involved  Green’s  functions.  The  ill-posedness  of 
the  problem  is  accentuated  due  to  the  involvement  of  both  G(r,rs)  and  G(rrf,r)  in  optical  imaging. 

Recently,  independent  component  analysis  (ICA)  has  been  shown  to  be  an  effective  method  to  recover 
unobserved  signals  (sources)  from  several  observed  mixtures. [5,  6]  In  this  report,  we  propose  a  novel  approach 
for  opt.ical  imaging.  The  distortion  of  the  light  exitance  is  interpreted  as  a  m.ixtuTe  of  virtual  sources  where 
the  virtual  source  and  the  mixing  matrix  represent  the  Green’s  functions  G(r,r5)  and  G(rrf,r)  respectively. 
The  mixing  matrix  and  the  virtual  sources  can  be  recovered  by  ICA  and  used  to  characterize  the  embedded 
inhomogeneities  by  fitting  the  individual  Green’s  functions. 

Foreign  absorptive  objects  located  inside  an  otherwise  uniform  infinite  slab  can  be  probed  using  CW  (or  time- 
resolved  or  frequency  modulated)  scanning  point  sources.  The  foreign  absorptive  objects  have  an  excessive 
absorption  (5//a(r)-  Denote  the  embedded  absorption  inhomogeneities  as  -  6^a{Li)cVi  where  c  is  the  speed 
of  light  in  the  medium  and  Vi  is  the  volume  of  the  ith  absorption  inhomogeneity,  the  scattered  wave  field 
can  be  written  as: 

Tl 

-</'sca{rd;rs)  =  y^G(rd,rj)gjG(rj,rs)  (1) 

where  G(r,r')  is  the  Green’s  function  of  the  uniform  background.  The  dependence  on  ui  is  not  explicitly 
shown  as  one  argument  of  the  Green’s  function  G  for  clarity. 

The  scattered  wave  (^sca  (1)  hence  represents  a  mixture  of  n  independent  virtual  sources  Sj{rs)  = 
qjG{Tj^rs)  by  the  mixing  matrix  A  whose  dj  entry  is  aj{Td)  =  G(rrf,rj),  and  reduced  to  the  well-known 
instantaneous  mixture  model: 

n 

x{r^)  =  As{t,)  =  (2) 

j=i 

where 


x(r,) 


(  0sca(r(ii !  Tj).  •  •  •  I 


T 


(3) 


s(r5)  =  (si(rs),...,Sn(rs))^ 

rp 

Bj  =  )>  •  •  •  j 

A  =  (a.1 , .  •  .  >  5 

and  rn  >  n  is  the  number  of  sensors.  Note  the  symmetry  about  the  source  and  the  sensors  in  Eq.  (1).  The 
role  of  sensors  and  sources  can  be  interchanged  freely  (the  reciprocal  property  of  light  propagation).  By 
using  the  mutual  independence  among  the  virtual  sources  Sj,  the  virtual  sources  sj  can  be  separated  up  to 
permutation  and  scaling.  The  virtual  sources  can  be  extracted  simultaneously  for  all  (a  parallel  scheme)  or 
one  by  one  (a  deflation  scheme). [5,  7] 

After  having  extracted  virtual  sources  and  their  corresponding  mixing  vectors  both  the  location  and 
strength  of  the  jth  object  can  be  computed  by  a  simple  fitting  procedure.  The  virtual  source  and  the 
mixing  vector  aj  are  the  scaled  version  of  the  Green’s  functions  G{rj,rs)  and  G{rd,Tj)  for  an  absorptive 
object  respectively,  i.e., 

Sj  =  ajG{rj,rs)  (4) 

SLj  =  PjGiTd^Vj), 

The  location  of  the  jth  object  and  the  strength  qj  of  the  jth  object  are  then  obtained  by  fitting  using 
Eq.  (4).  The  .strength  of  absorption  is  given  by  qj  —  ajPj, 

Fig.  (1)  displays  the  separation  and  characterization  of  two  simulated  point  absorbers,  each  of  an  absorption 
strength  of  unity,  inside  a  slab  of  thickness  50mm.  The  transport  mean  free  path  is  1mm.  The  incident  CW 
point  source  scans  throughout  a  set  of  21  x  21  grid  points  with  a  spacing  of  4.5mm  betwen  two  consecutive 
grid  points  on  one  side  of  the  slab.  The  light  exitance  on  the  other  side  of  the  slab  is  recorded  by  a  CCD 
camera  on  twice  denser  grid  points.  We  used  20%  additive  Gaussian  noise  in  this  simulation.  By  fitting  to 
Eq.  (4),  the  position  of  the  point  absorbers  are  found  to  be  very  close  to  the  input  values  with  an  en'or  less 
than  2mm  in  all  three  directions.  The  strengths  of  the  two  unit  absorbers  are  resolved  to  be  0.98  and  1.01 
with  an  error  of  less  than  2%,  respectively. 

Fig.  (2)  shows  the  result  of  unmixing  of  experimental  data  taken  for  two  parallel  absorptive  rods  immersed 
inside  an  lntralipid-10%  liquid  diluted  to  have  a  transport  mean  free  path  of  1mm.  The  thickness  of  the  slab 
is  50mm  and  the  scanning  step  is  5mm.  The  locations  are  found  to  be  a:  ==  —17.8mm  and  z  —  26.4mm  for 
the  first  rod,  and  x  ^  14.1mm  and  ^  ^  22.7ram  for  the  second  rod,  approximately  agreeing  with  the  input 
values. 

This  proposed  approach  is  able  to  locate  and  reconstruct  the  optical  properties  of  point  inhomogeneities, 
well  suited  for  the  cases  of  sparse  inhomogeneities.  In  practice,  the  inhomogeneity  of  a  finite  size  can  be 
approximated  by  a  point  one  as  long  as  its  dimension  is  much  smaller  than  the  distance  from  the  object  to 
the  surface  where  the  light  sources  and  detectors  are  placed.  The  same  approach  can  be  applied  to  mixed 
(absorptive  and  diffusing)  inhomogeneities.  The  shape  of  the  inhomogeneity  may  also  be  estimated  by  using 
a  similar  technique  to  back-propagation. 
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Fig.  1.  (a)  The  setup,  (b)  The  virtual  source  {G{r,  r^)),  its  mixing  vector  (G(rd,r))  and  the  fitting  to  Eq.  (4). 
The  left  column  is  for  the  first  point  source  at  (50,  60,30)mm  and  the  right  column  is  for  the  second  point 
source  at  (30, 30, 20)mm. 


Fig.  2.  (a)  The  setup,  (b)  The  virtual  sources  (^(r.r^))  and  the  mixing  vectors  (G(rd,r))  for  the  first  rod 
(the  upper  row)  and  the  second  rod  (the  lower  row). 
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Multiple  passages  of  light  through  an  absorption  inhomogeneity  of  finite  size  deep  within  a  turbid  medium 
are  analyzed  for  optical  imaging  by  use  of  the  self-energy  diagram.  The  nonlinear  correction  becomes  more 
important  for  an  inhomogeneity  of  a  larger  size  and  with  greater  contrast  in  absorption  with  respect  to  the  host 
background.  The  nonlinear  correction  factor  agrees  well  with  that  from  Monte  Carlo  simulations  for  cw  light. 
The  correction  is  approximately  50% -75%  in  the  near  infrared  for  an  absorption  inhomogeneity  with  the 
typical  optical  properties  found  in  tissues  and  five  times  the  size  of  the  transport  mean  free  path.  ©  2004 
Optical  Society  of  America 
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The  main  objective  of  optical  imaging  of  turbid  media 
is  to  locate  and  identify  the  embedded  inhomogeneities 
by  essentially  inverting  the  difference  in  photon 
transmittance  in  the  time  or  frequency  domains  due 
to  the  presence  of  these  inhomogeneities.^"^  The  key 
quantity  involved  is  the  Jacobian,  which  quantifies 
the  influence  on  the  detected  signal  due  to  the  change 
of  the  optical  parameters  of  the  medium.  The  linear 
perturbation  approach  is  suitable  for  calculating  the 
Jacobian  for  only  a  small  and  weak  absorption  in¬ 
homogeneity  and  is  not  valid  when  the  absorption 
strength  is  large.^  This  failure  can  be  attributed  to 
the  multiple  passages  through  the  abnormal  site  by 
the  photon.  The  most  important  correction  is  the 
self-energy  correction,®  which  takes  into  account  the 
repeated  visits  made  by  a  photon  through  the  site  up 
to  an  infinite  number  of  times.  The  presence  of  other 
inhomogeneity  islands  can  be  ignored  because  the 
photon  propagator  decreases  rapidly  with  the  distance 
between  two  separate  sites. 

In  this  Letter  the  nonlinear  correction  for  an  absorp¬ 
tion  inhomogeneity  of  a  large  strength  due  to  repeated 
visits  by  the  photon  is  modeled  by  a  nonlinear  correc¬ 
tion  factor  (NCF)  to  the  linear  perturbation  approach. 
The  NCF  as  a  function  of  the  size  and  the  strength 
of  the  inhomogeneity  is  estimated  by  use  of  the  self¬ 
energy  diagram.  The  NCF  is  obtained  from  the  cu- 
mulant  approximation  to  the  radiative  transfer  and 
verified  by  Monte  Carlo  simulations  for  cw  light.  The 
magnitude  of  the  NCF  is  0.5-1  for  an  absorptive  in¬ 
homogeneity  of  up  to  5lt  (It  is  the  mean  transport  free 
path  of  light)  and  of  the  typical  optical  properties  of 
human  tissues  (piah/c  ^  0.01-0.05,  where  pa  is  the 
absorption  coefficient  and  c  is  the  speed  of  light  in  the 
medium). 

If  we  consider  an  absorption  site  centered  at  r  and 
far  away  from  both  the  source  and  the  detector,  the 
change  in  the  detected  light  A7  at  the  detector  from 
a  modulated  point  source  at  including  the  multiple 
passages  through  the  site  is  given  by 


A/  =  -G(rd,ai|r)V5/x„(F)X[-Vseir(^^;R)V5M„(r)r 

n-0 

X  G(r,o;  Ir^) 

=  -G(rrf,  1  r)  7—^— — - — =7^77 — -- 

xG(r,a>|r,),  (1) 

where  Spa  is  the  excess  absorption  of  the  absorption 
site  of  size  R  and  volume  V,  oj  is  the  modulation  fre¬ 
quency  of  light,  G  is  the  propagator  of  photon  migra¬ 
tion  in  the  background  medium,  and 

Ns2\{(<o;R)=  ^  J^G(r2,a;  |ri)d^r2d'Vi  (2) 

is  the  self-propagator  that  describes  the  probability 
that  a  photon  revisits  volume  V,  Here  G(r2,  a>|ri) 
gives  the  probability  density  that  a  photon  leaves  the 
volume  at  rj  and  re-enters  it  at  r2.  The  scattering 
property  of  the  site  is  the  same  as  that  of  the  back¬ 
ground.  In  Eq.  ( 1)  G  (rd ,  I  r)  and  G  (F ,  |  rj)  are  well 
modeled  by  the  center-moved  diffusion  model  as  long 
as  separations  Irj  -  F |  and  Ir^  -  F|  are  much  greater 
than  Z/.’  However,  the  diffusion  Greenes  function  can¬ 
not  be  used  in  Eq.  (2)  to  evaluate  iVeeifCw;/?)  because 
the  diffusion  approximation  breaks  down  when  ri  is  in 
the  proximity  of  r2. 

Comparing  Eq.  (1)  with  the  standard  linear  pertur¬ 
bation  approach,  the  nonlinear  multiple  passage  effect 
of  an  absorption  site  is  represented  by  a  NCF: 

NCF  =  [1  +  N,,u((o;R)V8pa{r)r\  (3) 

This  factor  serves  as  a  universal  measure  of  the  non¬ 
linear  multiple-passage  effect  as  long  as  the  absorption 
site  is  far  from  both  the  source  and  the  detector  and  its 
size  is  much  smaller  than  its  distance  to  both  the  source 
and  the  detector.  This  correction  is  more  significant 
when  the  NCF  is  further  away  from  unity. 
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Photon  propagator  A7’(r2,  i  lri,s),  the  probability 
that  a  photon  propagates  from  position  ri  with  propa¬ 
gation  direction  s  to  position  r2  in  time  t,  for  any  sepa¬ 
ration  between  ri  and  r2,  was  recently  derived^’®  in  a 
form  of  the  cumulant  approximation  to  the  radiative 
transfer. 

In  the  case  of  interest  in  which  the  absorption  site 
is  deep  inside  the  medium,  the  photon  distribution 
is  isotropic.  The  photon  propagator  is  simplified 
to  Neff(r,^)  =  iVeff(lr2  -  ril.O,  which  is  obtained  by 
averaging  N(r 2,  H  ri ,  s)  over  the  propagation  direction 
s  of  light  over  the  in  solid  angle.  In  the  frequency 
domain  this  effective  propagator  is  approximately 
given  by 


NeairyO)) 


exp  I 


Airr'^c 

^  exp(-K/f) 

A'TvDrKlt 


("I””'"') 

sinh(/cr) , 


exp(-/cr) 

ATrDrKlt 


sinh(Klt) , 


r  <  //  >  (4) 


r>lt 


where  D  =  ltc/3  and  k  =  [3(/Xa  ~  whose 

sign  is  chosen  with  a  nonnegative  real  part.  The  two 
terms  in  Nctr  when  r  <  It  represent  ballistic  and  dif¬ 
fusion  contributions,  respectively.  The  ballistic  term 
does  not  depend  on  scattering  because  the  photon  dis¬ 
tribution  involved  is  already  isotropic.  Only  diffusion 
contributes  to  ATeff  when  r  >  Z^.  The  self-propagator 
for  an  absorption  sphere  deep  inside  the  medium  is 
given  by 


NsMi(o;R)  =  ^  -  ril,  w)d^r2d^ri 


1 

=  77  Kcr(r,(o)yo(r)47rr^dr,  (5) 

V  Jo 

where  7o(r)  =  1  -  (3r/4R)  +  {1/13)  (r/R)^  is  the 
characteristic  function  for  a  uniform  sphere.®  An  ab¬ 
sorption  site  of  an  arbitrary  shape  can  be  treated  the 
same  way.  The  exact  self-propagator  must  be  com¬ 
puted  by  a_ numerical  quadrature.  A  good  approxi¬ 
mation  of  Nse\f{(*>’yR)  is 

N,euico-.R)  =  jt 


+  +  f£l/2 

('i.  ^2  +  1  _  A  ^-1  +  _L  >  (6) 

V  5  ^  2  16  ^  320  ^  / 

+ o{k^),  i>m 


by  use  of  relation  (4),  where  f  ^  R/lt  when  \k\R  <$C  1. 
The  exact  and  approximate  versions  of  dimensionless 
self-propagator  A/’geifVZ^'^c  when  /<  =  0  are  plotted  as 
solid  and  dashed  curves,  respectively,  in  Fig.  1(a).  Di¬ 
mensionless  self-propagator  NseifVlJ'^c  depends  solely 
on  two  dimensionless  quantities  Klf  of  the  background 
and  R/lt  of  the  absorbing  sphere. 

It  is  worthwhile  to  point  out  that  self-propagator 
in  time  Nse\f{lyH),  tbe  inverse  Fourier  transform  of 


Eq.  (5),  includes  the  contribution  from  the  ballistic  mo¬ 
tion  of  the  photon  when  the  photon  passes  through  the 
site.  This  ballistic  contribution  manifests  itself  as  the 
linear  decay  of  Nse]dt;R)V  in  the  form  of  ro(cZ)  near 
the  origin  of  the  time,  followed  by  a  transition  to  dif¬ 
fusion  [Fig.  1(b)]. 

The  NCF  is  obtained  by  plugging  Eq.  (5)  or  (6)  into 
Eq.  (3).  In  particular,  we  have 


NCF- 


^  <1/2 


»  (7) 


^  >1/2 


where  q  =  V5/ia(r)/Zfc  is  the  dimensionless  strength 
of  the  absorber  when  \k\R  «  1.  For  an  absorber  of 
fixed  qr  >  0,  the  effectiveness  of  absorbing  light  is  di¬ 
minished  (the  NCF  decreases)  when  its  size  is  reduced. 
This  can  be  understood  from  the  fact  that  the  photon 
spends  less  time  per  volume  inside  the  absorber  of  a 
smaller  dimension  because  of  the  ballistic  motion  of  the 
photon  after  each  scattering  event.  The  photon  leaves 
a  small  site  {R  <  It)  in  an  almost  straight  line.  The 
diffusion  behavior  for  an  individual  photon  is  observed 
only  after  a  large  number  of  scattering  and  on  a  scale 
larger  than  Z;. 

Figure  2  shows  plots  of  the  NCF  versus  absorber 
size  for  typical  absorbers  of  excess  absorption  5  fMah/c 
equal  to  0.01  and  0.05.  The  nonlinear  correction  fac¬ 
tor  generally  decreases  with  the  size  of  the  absorber 
whose  excess  absorption  is  fixed.  With  the  increase 
of  the  background  absorption  and  the  modulation  fre¬ 
quency,  the  nonlinear  correction  becomes  less  accentu¬ 
ated.  The  phase  delay  is  larger  for  higher  modulation 
frequencies  and  less  background  absorption. 

Monte  Carlo  simulations^®  are  performed  for  cw  light 
propagating  in  a  uniform  nonabsorbing  and  isotropic 
scattering  slab.  The  thickness  of  the  slab  is  L  —  SOZ; . 
A  spherical  absorber  of  radius  R  is  located  at  the  center 
(0, 0,  L/2)  of  the  slab.  The  excess  absorption  of  the  ab¬ 
sorber  is  dfiaU/c  =  0.01.  The  absorber  has  the  same 
scattering  property  as  the  backgroimd.  The  details  of 
the  Monte  Carlo  computation  were  provided  in  a  pre¬ 
vious  publication.*^  The  correlated  sampling  method 
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(b) 


Fig.  1.  (a)  Self-propagator  Arseif(£t;;i2)VZf  *c  and  its  ap¬ 
proximation  form  when  /c  —  0.  (b)  Self-propagator  for 
spheres  of  various  radii  in  the  time  domain  inside  a 
nonabsorbing  medium. 
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Fig.  2.  NCF  (magnitude  and  phase  angle)  versus  the  size 
of  absorbers  whose  excess  absorption  8/iah/c  equals  0.01 
and  0.05.  Note  that  ^  3(/Za  -  iu))lt/c  for  the  back¬ 
ground  medium. 


Fig.  3.  (a)  Theoretical  nonlinear  correction  factors  from 

numerical  quadrature  (Exact),  the  approximate  form  of 
relation  (7)  (Approx),  and  Monte  Carlo  simulations  (MC). 
Results  from  four  independent  Monte  Carlo  simulations 
are  shown  for  each  radius.  The  standard  linear  perturba¬ 
tion  approach  corresponds  to  horizontal  line  NCF  =  1  (not 
shown  in  the  figure),  (b)  Percentage  change  of  the  cw 
transmittance  from  the  experimental  data  given  in  Fig.  9 
of  Ref.  5  compared  with  the  theoretical  predictions  made 
by  the  standard  linear  perturbation  approach  (StdPert) 
and  those  including  NCF  (Exact  and  Approx). 


is  used  in  each  simulation  to  reduce  variance. A 
single  simulation  is  used  to  compute  the  change  in  light 
transmittance  due  to  the  presence  of  the  absorption 
site  and  the  corresponding  NCF.  Figure  3(a)  shows 
the  NCFs  obtained  from  numerical  quadrature,  the  ap¬ 
proximate  form  of  relation  (7),  and  Monte  Carlo  simu¬ 
lations.  The  agreement  between  our  theoretical  NCF 
and  that  from  Monte  Carlo  simulations  is  excellent. 
The  slight  difference  between  them  at  large  radii  is 
accounted  for  by  the  fact  that  the  sphere  can  no  longer 
be  regarded  as  small  compared  with  the  dimensions 
of  the  slab.  The  probability  of  a  photon  revisiting  a 
large  sphere  is  overestimated  by  Eq.  (5)  for  the  sphere 
located  at  the  center  of  the  slab.^^ 


Figure  3(b)  shows  the  percentage  change  of  the  cw 
transmittance  estimated  from  the  experimental  data 
given  in  Fig.  9  of  Ref.  5.  The  relevant  parameters 
of  the  experiment  are  summarized  in  the  inset.  The 
theoretical  predictions  from  the  linear  perturbation  ap¬ 
proach  with  and  without  the  nonlinear  correction  are 
also  shown  in  Fig.  3(b),  assuming  a  collimated  point 
source  and  a  point  detector  in  a  confocal  setup.  Our 
theoretical  prediction  with  nonlinear  correction  pro¬ 
vides  a  significant  improvement  over  linear  pertur¬ 
bation  and  agrees  much  better  with  the  experimental 
result. 

The  typical  value  of  the  absorption  coefficient  of 
human  tissues  in  the  near  infrared  indicates  that 
f^ahlc  0.01-0.05. This  fact  should  put  our 
results  on  NCFs  in  this  range  (Figs.  2  and  3)  into 
perspective.  The  nonlinear  correction  becomes  more 
important  for  an  inhomogeneity  of  a  larger  size  and 
with  greater  contrast  in  absorption  with  respect  to 
the  background.  The  value  of  the  NCF  decreases 
from  --'0.75  to  —0.5  for  an  absorption  site  of  radius  5lt 
with  excess  absorption  Sfiah/c  increasing  from  0.01  to 
0.05.  The  standard  linear  perturbation  approach  in 
optical  imaging  should  be  augmented  to  include  this 
nonlinear  correction. 
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